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Abstract 

The current status of numerical solutions for the equations of ideal general rel- 
ativistic hydrodynamics is reviewed. With respect to an earlier version of the 
article the present update provides additional information on numerical schemes 
and extends the discussion of astrophysical simulations in general relativistic hy- 
drodynamics. Different formulations of the equations are presented, with special 
mention of conservative and hyperbolic formulations well-adapted to advanced 
numerical methods. A large sample of available numerical schemes is discussed, 
paying particular attention to solution procedures based on schemes exploiting 
the characteristic structure of the equations through linearized Riemann solvers. 
A comprehensive summary of astrophysical simulations in strong gravitational 
fields is presented. These include gravitational collapse, accretion onto black 
holes and hydrodynamical evolutions of neutron stars. The material contained 
in these sections highlights the numerical challenges of various representative 
simulations. It also follows, to some extent, the chronological development of 
the field, concerning advances on the formulation of the gravitational field and 
hydrodynamic equations and the numerical methodology designed to solve them. 



1 Introduction 

The description of important areas of modern astronomy, such as high-energy 
astrophysics or gravitational wave astronomy, requires General Relativity. High 
energy radiation is often emitted by highly relativistic events in regions of strong 
gravitational fields near compact objects such as neutron stars or black holes. 
The production of relativistic radio jets in active galactic nuclei, explained by 
pure hydrodynamical effects as in the twin-exhaust model [38], by hydromag- 
netic centrifugal acceleration as in the Blandford-Payne mechanism [37], or by 
electromagnetic extraction of energy as in the Blandford-Znajek mechanism [39], 
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involves an accretion disk around a rotating supermassive black hole. The dis- 
covery of kHz quasi-periodic oscillations in low-mass X-ray binaries extended 
the frequency range over which these oscillations occur into timescales associ- 
ated with the relativistic, innermost regions of accretion disks (see, e.g. [289]). A 
relativistic description is also necessary in scenarios involving explosive collapse 
of very massive stars SOA^q) to a black hole (in the so-called coUapsar and 
hypernova models), or during the last phases of the coalescence of neutron star 
binaries. These catastrophic events are believed to exist at the central engine 
of highly energetic 7-ray bursts (GRBs) [216, 202, 307, 217]. In addition, non 
spherical gravitational collapse leading to black hole formation or to a supernova 
explosion, and neutron star binary coalescence are among the most promising 
sources of detectable gravitational radiation. Such astrophysical scenarios con- 
stitute one of the main targets for the new generation of ground-based laser 
interferometers, just starting their gravitational wave search (LIGO, VIRGO, 
GEO600, TAMA) [286, 207]. 

A powerful way to improve our understanding of the above scenarios is 
through accurate, large scale, three-dimensional numerical simulations. Nowa- 
days, computational general relativistic astrophysics is an increasingly impor- 
tant field of research. In addition to the large amount of observational data 
by high-energy X- and 7-ray satellites such as Chandra, XMM-Newton or IN- 
TEGRAL, and the new generation of gravitational wave detectors, the rapid 
increase in computing power through parallel supercomputers and the associ- 
ated advance in software technologies is making possible large scale numerical 
simulations in the framework of general relativity. However, the computational 
astrophysicist and the numerical relativist face a daunting task. In the most 
general case, the equations governing the dynamics of relativistic astrophysical 
systems are an intricate, coupled system of time-dependent partial differen- 
tial equations, comprising the (general) relativistic (magneto-) hydrodynamic 
(MHD) equations and the Einstein gravitational field equations. In many cases, 
the number of equations must be aTigmcmtcd to account for non-adiabatic pro- 
cesses, e.g., radiative transfer or sophisticated microphysics (realistic equations 
of state for nuclear matter, nuclear physics, magnetic fields, etc.). 

Nevertheless, in some astrophysical situations of interest, e.g., accretion of 
matter onto compact objects or oscillations of relativistic stars, the 'test-fluid' 
approximation is enough to get an accurate description of the underlying dy- 
namics. In this approximation the fluid self-gravity is neglected in comparison 
to the background gravitational field. This is best exemplified in accretion prob- 
lems where the mass of the accreting fiuid is usually much smaller than the 
mass of the compact object. Additionally, a description employing ideal hydro- 
dynamics (i.e., with the stress-energy tensor being that of a perfect fluid), is 
also a fairly standard choice in numerical astrophysics. 

The main purpose of this review is to summarize the existing efforts to solve 
numerically the equations of (ideal) general relativistic hydrodynamics. To this 
aim, the most important numerical schemes will be first presented in some de- 
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tail. Prominence will be given to the so-called Godunov-type schemes written in 
conservation form. Since [166] it has been gradually demonstrated ([97, 82, 244, 
87, 25, 297, 230]) that conservative methods exploiting the hyperbolic character 
of the relativistic hydrodynamic equations are optimally suited for accurate nu- 
merical integrations, even well inside the ultrarelativistic regime. The explicit 
knowledge of the characteristic speeds (eigenvalues) of the equations, together 
with the corresponding eigenvectors, provides the mathematical (and physical) 
framework for such integrations, by means of either exact or approximate Rie- 
mann solvers. 

The article includes, furthermore, a comprehensive description of 'relevant' 
numerical applications in relativistic astrophysics, including gravitational col- 
lapse, accretion onto compact objects and hydrodynamical evolution of neu- 
tron stars. Numerical simulations of strong- field scenarios employing Newtonian 
gravity and hydrodynamics, as well as possible post-Newtonian extensions, have 
received considerable attention in the literature and will not be covered in the 
review, which focuses in relativistic simulations. Nevertheless, we must empha- 
size that most of what is known about hydrodynamics near compact objects, 
in particular in black hole astrophysics, has been accurately described using 
Newtonian models. Probably the best known example is the use of a pseudo- 
Newtonian potential for non-rotating black holes which mimics the existence 
of an event horizon at the Schwarzschild gravitational radius [218], which has 
allowed accurate interpretations of observational phenomena. 

The organization of the article is as follows: Section 2 presents the equations 
of general relativistic hydrodynamics, summarizing the most relevant theoretical 
formulations which, to some extent, have helped to drive the development of nu- 
merical algorithms for their solution. Section 3 is mainly devoted to describing 
numerical schemes specifically designed to solve non-linear hyperbolic systems 
of conservation laws. Hence, particular emphasis will be paid on conservative 
high-resolution shock-capturing (HRSC) upwind methods based on linearized 
Ricmann solvers. Alternative schemes such as Smoothed Particle Hydrodynam- 
ics (SPH), (pseudo-) spectral methods and others will be briefly discussed as 
well. Section 4 summarizes a comprehensive sample of hydrodynamical sim- 
ulations in strong-field genca'al relativistic astrophysics. Finally, in Section 5 
we provide additional technical information needed to build up upwind HRSC 
schemes for the general relativistic hydrodynamics equations. Geometrized units 
(G = c = 1) are used throughout the paper except where explicitly indicated, 
as well as the metric conventions of [187]. Greek (Latin) indices run from to 
3 (1 to 3). 



3 



2 Equations of general relativistic hydrodynam- 
ics 



The general relativistic hydro dynamic equations consist of the local conservation 
laws of the stress-energy tensor, T'^" (the Bianchi identities) and of the matter 
current density, J** (the continuity equation) 

V^T'*'' = 0, (1) 

V^J^ = 0. (2) 

As usual stands for the covariant derivative associated with the four- 
dimensional spacctimc metric g^y. The density current is given by = pu^, 

representing the fluid 4-velocity and p the rest-mass density in a locally 
inertial reference frame. 

The stress-energy tensor for a non-perfect fluid is defined as 

T^"' = p(l + e)uV + {p- C0)h^"' - 2r](7^"' + qV + q-'u" (3) 

where e is the rest frame specific internal energy density of the fluid, p is the 

pressure and ft,^" is the spatial projection tensor /i^"^ = u^v^ + g^^^ . In addition, 
•q and C, are the shear and bulk viscosities. The expansion 9, describing the 
divergence or convergence of the fluid world lines is defined a,s 6 = V^m'*. The 
symmetric, trace- free, and spatial shear tensor a'^" , is defined by 

a^"" = hv^u^h""" + VocU^h"'^') - \eh'"', (4) 

and, finally, g** is the energy fiux vector. 

In the following we will neglect non-adiabatic effects, such as viscosity or 
heat transfer, assuming the stress-energy tensor to be that of a perfect fiuid 

T'"' = phuV + pg^"", (5) 

where we have introduced the relativistic specific enthalpy, h, defined by 

h = l + e+-. (6) 

P 

Introducing an explicit coordinate chart (a;°,x*) the previous conservation 
equations read 



d 



^T"- = -^T';^T'^\ (8) 



where the scalar x° represents a foliation of the spacetime with hypersurfaces 
(coordinatized by a;'). Additionally, is the volume element associated with 
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the 4-metric, with g = det(5jui/), and FJ^^ are the 4-dimensional Christoffel 
symbols. 

In order to close the system, the equations of motion (1) and the continuity 
equation (2) must be supplemented with an equation of state (EOS) relating 
some fundamental thermodynamical quantities. In general, the EOS takes the 
form 

P = p{p.e)- (9) 

Due to their simplicity, the most widely employed EOS in numerical simu- 
lations are the ideal fluid EOS, p = {T — l)ps, where F is the adiabatic index, 
and the polytropic EOS (e.g. to build equilibrium stellar models), p = Kp^ = 
^pi+i/jv^ being the polytropic constant and N the polytropic index. 

In the 'test-fluid' approximation, where the fluid self-gravity is neglected, the 
dynamics of the system is completely governed by Eqs. (1) and (2), together with 
the EOS (9). In those situations where such approximation does not hold, the 
previous equations must be solved in conjunction with the Einstein gravitational 
fleld equations, 

G^"' = SttT'"', (10) 

which describe the evolution of the geometry in a dynamical spacetime. A 
detailed description of the various numerical approaches to solve the Einstein 
equations is beyond the scope of the present article (see, e.g. Lehncr [154] for a 
recent review) . We only mention that the formulation of the Einstein equations 
as an initial value (Cauchy) problem, in the presence of matter fields, adopting 
the so-called 3-1-1 decomposition of the spacetime [19] can be found in, e.g., [315]. 
Given a choice of gauge, the Einstein equations in the 3-1-1 formalism [19] split 
into evolution equations for the 3-metric jij and the extrinsic curvature Kij 
(the second fundamental form), and constraint equations, the Hamiltonian and 
momentum constraints, that must be satisfied at every time slice. Long-term 
stable evolutions of the Einstein equations have recently been accomplished us- 
ing various reformulations of the original 3-1-1 system (see, e.g. [29, 257, 8, 93] for 
simulations involving matter sources, and [11] and references therein for (vac- 
uum) black hole evolutions). Alternatively, a characteristic initial value problem 
formulation of the Einstein equations was developed in the 1960s by Bondi, van 
der Burg and Metzner [48], and Sachs [247]. This approach has gradually ad- 
vanced to a state where long-term stable evolutions of caustic-free spacctimes 
in multidimcnsions are possible, even including matter fields (see [154] and ref- 
erences therein). A recent review of the characteristic formulation is presented 
in a Living Reviews article by Winicour [305]. Examples of this formulation 
in general relativistic hydrodynamics are discussed in various sections of the 
present article. 

Traditionally, most of the approaches for numerical integrations of the gen- 
eral relativistic hydrodynamic equations have adopted spacelike foliations of the 
spacetime, within the 3-1-1 formulation. Recently, however, covariant forms of 
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these equations, well suited for advanced numerical methods, have also been 
developed. This is reviewed next in a chronological way. 



2.1 Spacelike (3+1) approaches 

In the 3+1 (ADM) formulation [19], spacetime is foliated into a set of non- 
intersecting spacelike hypersurfaces. There are two kinematic variables describ- 
ing the evolution between these surfaces: the lapse function a. which describes 
the rate of advance of time along a timelike unit vector normal to a surface, 
and the spacelike shift vector /3* that describes the motion of coordinates within 
a surface. 

The line element is written as 



where jij is the 3-metric induced on each spacelike slice. 

2.1.1 1+1 Lagrangian formulation (May and White) 

The pioneering numerical work in general rclativistic hydrodynamics dates back 
to the one-dimensional gravitational collapse code of May and White [173, 174]. 
Building on theoretical work by Misner and Sharp [186], May and White devel- 
oped a time-dependent numerical code to solve the evolution equations describ- 
ing adiabatic spherical collapse in general relativity. This code was based on a 
Lagrangian finite difference scheme (see Section 3.1), in which the coordinates 
are co-moving with the fluid. Artificial viscosity terms were included in the 
equations to damp the spurious numerical oscillations caused by the presence 
of shock waves in the flow solution. May and White's formulation became the 
starting point of a large number of numerical investigations in subsequent years 
and, hence, it is worth describing its main features in some detail. 

For a spherically symmetric spacetime the line element can be written as 



ds^ = -a2(m, t)dt'^ + b^{m, t)dm^ + R^{rn, t){dd^ + sin^ Odc/)^), (12) 



m being a radial (Lagrangian) coordinate, indicating the total rest-mass en- 
closed inside the circumference 2TrR{m,t). 

The co-moving character of the coordinates leads, for a perfect fluid, to a 
stress-energy tensor of the form: 



In these coordinates the local conservation equation for the baryonic mass, 
Eq. (2), can be easily integrated to yield the metric potential b: 



ds"^ = -(^2 - 3,B')dx"dx" + 2f3,dx'dx' 



+ ^ijdx^dx^ 



(11) 



Tl = T| = T| = -p, 
T^ = {l + e)p, 
T^ = if iii-v. 



(13) 
(14) 
(15) 



b = 



1 



(16) 



A-npB? ' 
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The gravitational field equations, Eq. (10), and the equations of motion, 
Eq. (1), reduce to the following quasi-linear system of partial differential equa- 
tions (see also [186]): 

dp 1 da , „ ,^„, 
^ + -—ph = Q, 17 
dm a am 



dt "^dt \p 
du ( ^_ dp 2r , M 



1 dpR- _ d^/drn^ ^^0) 



dt dR/dm' 



with the definitions u = a.nd T = l§^, satisfying = 1 - _ 2M 

Additionally, 

M= / 4nR'^p(l + e)—dm, (21) 
Jo dm 

represents the total mass interior to radius m at time t. The final system, 
Eqs. (17)-(20), is closed with an EOS of the form given by Eq. (9). 

Hydrodynamics codes based on the original formulation of May and White 
and on later versions (e.g., [294]) have been used in many non- linear simula- 
tions of supernova and neutron star collapse (see, e.g., [185, 280] and refer- 
ences therein), as well as in perturbativc computations of spherically symmetric 
gravitational collapse within the framework of the linearized Einstein equa- 
tions [251, 252]. In Section 4.1.1 below some of these simulations are discussed 
in detail. An interesting analysis of the above formulation in the context of 
gravitational collapse is provided by Miller and Sciama [183]. By comparing 
the Newtonian and relativistic equations, these authors showed that the net 
acceleration of the infalling mass shells is larger in general relativity than in 
Newtonian gravity. The Lagrangian character of May and White's formulation, 
together with other theoretical considerations concerning the particular coordi- 
nate gauge, has prevented its extension to multidimensional calculations. How- 
ever, for one-dimensional problems, the Lagrangian approach adopted by May 
and White has considerable advantages with respect to an Eulerian approach 
with spatially fixed coordinates, most notably the lack of numerical diffusion. 



2.1.2 3+1 Eulerian formulation (Wilson) 

The use of Eulerian coordinates in multidimensional numerical relativistic hy- 
drodynamics started with the pioneering work by Wilson [299] . Introducing the 
basic dynamical variables D, and E, representing the relativistic density, 
momenta and energy, respectively, defined as 

D = pu°, Si^ = phu^^vP, E = peu^, (22) 
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the equations of motion in Wilson's formulation [299, 300] are: 




(23) 




0, (24) 



(25) 



with the "transport velocity" given by = u^^ /u^ . We note that in the original 
formulation [300] the momentum density equation, Eq. (24), is only solved for 
the three spatial components, 5*^, and is obtained through the 4- velocity 
normalization condition u^u^ = — 1. 

A direct inspection of the system shows that the equations arc written as a 
coupled set of advection equations. In doing so, the terms containing derivatives 
(in space or time) of the pressure arc treated as source terms. This approach, 
hence, sidesteps an important guideline for the formulation of non-linear hyper- 
bolic systems of equations, namely the preservation of their conservation form. 
This is a necessary condition to guarantee correct evolution in regions of sharp 
entropy generation (i.e., shocks). Furthermore, some amount of numerical dissi- 
pation must be used to stabilize the solution across discontinuities. In this spirit 
the first attempt to solve the equations of general relativistic hydrodynamics in 
the original Wilson's scheme [299] used a combination of finite difference up- 
wind techniques with artificial viscosity terms. Such terms adapted the classic 
treatment of shock waves introduced by von Neumann and Richtmyer [296] to 
the relativistic regime (see Section 3.1.1). 

Wilson's formulation has been widely used in hydrodynamical codes de- 
veloped by a variety of research groups. Many different astrophysical sce- 
narios were first investigated with these codes, including axisymmetric stel- 
lar core-collapse [196, 194, 200, 26, 276, 229, 83], accretion onto compact ob- 
jects [125, 227], numerical cosmology [57, 58, 16] and, more recently, the coa- 
lescence of neutron star binaries [303, 304, 172]. This formalism has also been 
employed, in the special relativistic limit, in numerical studies of heavy-ion col- 
lisions [302, 176]. We note that in most of these investigations, the original 
formulation of the hydrodynamic equations was slightly modified by re-defining 
the dynamical variables, Eq. (22), with the addition of a multiplicative a factor 
(the lapse function) and the introduction of the Lorentz factor, W = avP: 



As mentioned before, the description of the evolution of self-gravitating mat- 
ter fields in general relativity requires a joint integration of the hydrodynamic 
equations and the gravitational field equations (the Einstein equations). Using 
Wilson's formulation for the fluid dynamics, such coupled simulations were first 



D = pW, 5^ = phWu, 



E = peW. 



(26) 
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considered in [300] , building on a vacuum, numerical relativity code specifically 
developed to investigate the head-on collision of two black holes [273] . The re- 
sulting code was axially symmetric and aimed to integrate the coupled set of 
equations in the context of stellar core collapse [86] . 

More recently, Wilson's formulation has been applied to the numerical study 
of the coalescence of binary neutron stars in general relativity [303, 304, 172] (see 
Section 4.4.2). These studies adopted an approximation scheme for the gravita- 
tional field, by imposing the simplifying condition that the three-geometry (the 
3-metric 7ij) is conformally flat. The line element, Eq. (11), then reads 

ds'^ = -(a^ - /3ip')dx°dx° + 2l3idx'dx° + (f)%jdx'dx\ (27) 

The curvature of the three metric is then described by a position dependent 
conformal factor times a flat-space Kronecker delta. Therefore, in this ap- 
proximation scheme all radiation degrees of freedom are removed, while the field 
equations reduce to a set of five Poisson-likc elliptic equations in flat spacetime 
for the lapse, the shift vector and the conformal factor. While in spherical sym- 
metry this approach is no longer an approximation, being identical to Einstein's 
theory, beyond spherical symmetry its quality degrades. In [142] it was shown 
by means of numerical simulations of extremely relativistic disks of dust that it 
has the same accuracy as the first post-Newtonian approximation. 

Wilson's formulation showed some limitations in handling situations involv- 
ing ultrarelativistic fiows {W » 2), as first pointed out by Centrella and Wil- 
son [58]. Norman and Winkler [209] performed a comprehensive numerical as- 
sessment of such formulation by means of special relativistic hydrodynamical 
simulations. Figure 1 reproduces a plot from [209] in which the relative error of 
the density compression ratio in the so-called relativistic shock reflection prob- 
lem - the heating of a cold gas which impacts at relativistic speeds with a solid 
wall and bounces back - is displayed as a function of the Lorentz factor W of 
the incoming gas. The source of the data is Ref. [58]. This flgure shows that 
for Lorentz factors of about 2 {v !v 0.86c), the threshold of the ultrarelativistic 
limit, the relative errors are between 5% and 7% (depending on the adiabatic 
exponent of the gas), showing a linear growth with W. 

Norman and Winkler [209] concluded that those large errors were mainly 
due to the way in which the artiflcial viscosity terms are included in the nu- 
merical scheme in Wilson's formulation. These terms, commonly called Q in 
the literature (see Section 3.1.1), are only added to the pressure terms in some 
cases, namely at the pressure gradient in the source of the momentum equation, 
Eq. (24), and at the divergence of the velocity in the source of the energy equa- 
tion. Eq. (25). However, [209] proposed to add the Q terms in a relativistically 
consistent way, in order to consider the artificial viscosity as a real viscosity. 
Hence, the hydrodynamic equations should be rewritten for a modified stress- 
energy tensor of the following form: 

T^"' = p{l+e+{p + Q)/p)uV + {p + Q)gi'r (28) 
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Figure 1: Results for the shock heating test of a cold, relativistically inflow- 
ing gas against a wall using the explicit Eulcrian techniques of Centrella and 
Wilson [58]. Dependence of the relative errors of the density compression ratio 
versus the Lorentz factor W for two different values of the adiabatic index of 
the flow, r = 4/3 (triangles) and F = 5/3 (circles) gases. The relative error is 
measured with respect to the average value of the density over a region in the 
shocked material. The data are from Centrella and Wilson [58] and the plot 
reproduces a similar one from Norman and Winkler [209] . 



In this way, for instance, the momentum equation takes the following form (in 
flat spacetime): 

^[{ph + Q)W'Vj] + ^[{ph + Q)W%V^] + ^-^±^ = 0. (29) 

In Wilson's original formulation Q is omitted in the two terms containing the 
quantity ph. In general Q is a non-linear function of the velocity and, hence, 
the quantity QW^V in the momentum density of Eq. (29) is a highly non- linear 
function of the velocity and its derivatives. This fact, together with the explicit 
presence of the Lorentz factor in the convective terms of the hydrodynamic 
equations, as well as the pressure in the specific enthalpy, make the relativistic 
equations much more coupled than their Newtonian counterparts. As a result 
Norman and Winkler proposed the use of implicit schemes as a way to describe 
more accurately such coupling. Their code, which in addition incorporates an 
adaptive grid, reproduces very accurate results even for ultra relativistic flows 
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with Lorentz factors of about 10 in one-dimensional, flat spacetime, simulations. 

Very recently, Anninos and Fragile [17] have compared state-of-the-art ar- 
tificial viscosity schemes and high-order non-oscillatory central schemes (see 
Section 3.1.3) using Wilson's formulation for the former class of schemes and 
a conservative formulation (similar to the one considered in [222, 220]; Section 
2.2.2) for the latter. These authors found, using a three-dimensional Cartesian 
code, that earlier results for artificial viscosity schemes in shock tube tests or 
shock reflection tests are not improved, i.e. the numerical solution becomes 
increasingly imstablc for shock velocities greater than about ^ 0.95c. On the 
other hand, results for the shock reflection problem with a second-order finite 
difference central scheme show the suitability of such a scheme to handle ul- 
trarelativistic flows, the underlying reason being, most likely, the use of a con- 
servative formulation of the hydrodynamic equations rather than the particular 
scheme employed (see Section 3.1.3). Tests concerning spherical accretion on 
to a Schwarzschild black hole using both schemes yield the maximum relative 
errors near the event horizon, as large as ~ 24% for the central scheme. 

2.1.3 3+1 conservative Eulerian formulation (Ibanez and coworkers) 

In 1991, Martf, Ibaficz and Mirallcs [166] presented a new formulation of the 
(Eulerian) general relativistic hydrodynamic equations. This formulation was 
aimed to take fundamental advantage of the hyperbolic and conservative charac- 
ter of the equations, contrary to the one discussed in the previous Section. From 
the numerical point of view, the hyperbolic and conservative nature of the rela- 
tivistic Euler equations allows for the use of schemes based on the characteristic 
fields of the system, translating to relativistic hydrodynamics existing tools of 
classical fluid dynamics. This procedure departs from earlier approaches, most 
notably in avoiding the need for artiflcial dissipation terms to handle discontin- 
uous solutions [299, 300] as well as implicit schemes as proposed in [209]. 

If a numerical scheme written in conservation form converges, it automati- 
cally guarantees the correct Rankine-Hugoniot (jump) conditions across discon- 
tinuities - the shock-capturing property (see, e.g. [155]). Writing the relativistic 
hydrodynamic equations as a system of conservation laws, identifying the suit- 
able vector of unknowns and building up an approximate Riemann solver per- 
mitted the extension of state-of-the-art high-resolution shock- capturing schemes 
(HRSC in the following) from classical fluid dynamics into the realm of relativ- 
ity [166]. 

Theoretical advances on the mathematical character of the relativistic hy- 
drodynamic equations were flrst achieved studying the special relativistic limit. 
In Minkowski spacetime, the hyperbolic character of relativistic hydrodynam- 
ics and magneto-hydrodynamics (MHD) was exhaustively studied by Anile and 
collaborators (see [14] and references therein) by applying Friedrichs' definition 
of hyperbolicity [104] to a quasi-linear form of the system of hydrodynamic 
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equations, 

^''(w)g^=0, (30) 

where A'^ are the Jacobian matrices of the system and w is a suitable set of 
primitive (physical) variables (see below). The system (30) is hyperbolic in the 
time direction defined by the vector field ^ with ^^^^ = —1, if the following 
two conditions hold: (i) det(^^^^) ^ and (ii) for any ( such that Cm^^ = 0, 
(ij,^^ = 1, the eigenvalue problem »4^(C/t — A^^)r = has only real eigenvalues 
{A„}„=i,...,5 and a complete set of right-eigenvectors {r„}„=i,...,5. Besides veri- 
fying the hyperbolic character of the relativistic hydrodynamic equations, Anile 
and collaborators [14] obtained the explicit expressions for the eigenvalues and 
eigenvectors in the local rest frame, characterized by = 6q. In Font et al. [97] 
those calciilations were extended to an arbitrary reference frame in which the 
motion of the fluid is described by the 4-velocity u'' = W{1, v^). 

The approach followed in [97] for the equations of special relativistic hy- 
drodynamics was extended to general relativity in [25] . The choice of evolved 
variables {conserved quantities) in the 3+1 Eulerian formulation developed by 
Banyuls et al. [25] differs slightly from that of Wilson's formulation [299] . It 
comprises the rest-mass density {D), the momentum density in the j-direction 
{Sj), and the total energy density (i5), measured by a family of observers which 
are the natural extension (for a generic spacetime) of the Eulerian observers 
in classical fluid dynamics. Interested readers are addressed to [25] for their 
definition and geometrical foundations. 

In terms of the so-called primitive variables w = {p,Vi,e), the conserved 
quantities are written as: 

D = pW, (31) 

Sj=phW^Vj, (32) 

E = phW'^ - p, (33) 

where the contravariant components = ^''vj of the three- velocity are defined 
as 

au^ a 

and W is the relativistic Lorentz factor W = au° = (1 — t;^)~^/^ with v'^ = 

With this choice of variables the equations can be written in conservation 
form. Strict conservation is only possible in flat spacetime. For curved space- 
times there exist source terms, arising from the spacetime geometry. However, 
these terms do not contain derivatives of stress-energy tensor components. More 
precisely, the first-order flux-conservative hyperbolic system, well suited for nu- 
merical applications, reads: 

4,(^.^^)-(w,. 



12 



with g = det{gixu) satisfying ^/—g = a^/j with 7 = det{jij). The state vector 
is given by 

Viw) = {D,Sj,T), (36) 
with T = E — D. The vector of fluxes is 

F^(w) = i^D (v^ ^) ,S, (v^ - ^) +pS),r(v^ - ^) + pv^^ , (37) 

and the corresponding sources S(w) are 

S(w) = (0,T- - rt^gs,^ , a (t-°|^ - r-r°,)) . (38) 

The local characteristic structure of the previous system of equations was 
presented in [25]. The eigenvalues (characteristic speeds) of the corresponding 
Jacobian matrices are all real (but not distinct, one showing a threefold degen- 
eracy as a result of the assumed directional splitting approach) and a complete 
set of right-eigenvectors exists. System (35) satisfies, hence, the definition of hy- 
perbolicity. As it will become apparent in Section 3.1.2 below, the knowledge of 
the spectral information is essential in order to construct HRSC schemes based 
on Riemann solvers. This information can be found in [25] (see also [100]). 

The range of applications considered so far in general relativity employing 
the above formulation of the hydrodynamic equations, Eq. (35)-(38), is still 
small and mostly devoted to the study of stellar core collapse and accretion flows 
on to black holes (see Sections 4.1.1 and 4.2 below). In the special relativistic 
limit this formulation is being successfully applied to simulate the evolution 
of (ultra-) relativistic extragalactic jets, using numerical models of increasing 
complexity (sec, e.g. [170, 12]). The first applications in general relativity were 
performed, in one spatial dimension, in [166], using a slightly different form of 
the equations. Preliminary investigations of gravitational stellar collapse were 
attempted by coupling the above formulation of the hydrodynamic equations 
to a hyperbolic formulation of the Einstein equations developed by [42] . These 
results arc discussed in [164, 41]. More recently, successful evolutions of fully 
dynamical spacetimcs in the context of adiabatic stellar core-collapse, both in 
spherical symmetry and in axisymmetry, have been achieved [132, 244, 71]. 
These investigations are considered in Section 4.1.1 below. 

An ambitious three-dimensional, Eulerian code which evolves the coupled 
system of Einstein and hydrodynamics equations was developed by Font et 
al. [100] (see Section 3.3.2). The formulation of the hydrodynamic equations in 
this code follows the conservative Eulerian approach discussed in this section. 
The code is constructed for a completely general spacetime metric based on a 
Cartesian coordinate system, with arbitrarily specifiable lapse and shift condi- 
tions. In [100] the spectral decomposition (eigenvalues and right-eigenvectors) 
of the general relativistic hydrodynamic equations, valid for general spatial met- 
rics, was derived, extending earlier results of [25] for non-diagonal metrics. A 
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complete set of left-eigenvectors was presented by Ibanez et al. [130]. Due to 
the paramount importance of the characteristic structure of the equations in the 
design of upwind HRSC schemes based upon Riemann solvers, we summarize 
all necessary information in Section 5.2 of the article. 

2.2 Covariant approaches 

General (covariant) conservative formulations of the general relativistic hydro- 
dynamic equations for ideal fluids, i.e., not restricted to spacelike foliations, 
have been reported in [82] and, more recently, in [222, 220]. The form invari- 
ance of these approaches with respect to the nature of the spacetime foliation 
implies that existing work on highly specialized techniques for fluid dynamics 
(i.e. HRSC schemes) can be adopted straightforwardly. In the next two sections 
we describe the existing covariant formulations in some detail. 

2.2.1 Eulderink and Mellema 

Eulderink and Mellema [80, 82] first derived a covariant formulation of the 
general relativistic hydrodynamic equations. As in the formulation discussed in 
Section 2.1.3, these authors took special care of the conservative form of the 
system, with no derivatives of the dependent fluid variables appearing in the 
source terms. Additionally, this formulation is strongly adapted to a particular 
numerical method based on a generalization of Roe's approximate Riemann 
solver. Such solver was first applied to the non-relativistic Euler equations in 
[242] and has been widely employed since in simulating compressible flows in 
computational fluid dynamics. Furthermore, their procedure is specialized for a 
perfect fluid EOS, p = {T — l)pe, T being the (constant) adiabatic index of the 
fluid. 

After the appropriate choice of the state vector variables, the conservation 
laws, Eqs. (7) and (8), are re- written in flux-conservative form. The flow vari- 
ables are then expressed in terms of a parameter vector co as 



where = Ku", = K-^ and K'^ = ^/^ph = —gapto^LO^. The vec- 
tor F° represents the state vector (the unknowns), and each vector F* is the 
corresponding flux in the coordinate direction x*. 

Eulderink and Mellema computed the exact "Roe matrix" [242] for the vec- 
tor (39) and obtained the corresponding spectral decomposition. The character- 
istic information is used to solve the system numerically using Roe's generalized 
approximate Riemann solver. Roe's linearization can be expressed in terms of 
the average state w = ^^^^^ , where L and R denote the left and right states in 
a Riemann problem (see Section 3.1.2). Further technical details can be found 
in [80, 82]. 



-pa _ 




) 



(39) 
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The performance of this general relativistic Roe solver was tested in a number 
of one-dimensional problems for which exact solutions are known, including non- 
relativistic shock tubes, special relativistic shock tubes and spherical accretion 
of dust and a perfect fluid onto a (static) Schwarzschild black hole. In its special 
relativistic version it has been used in the study of the confinement properties 
of relativistic jets [81]. However, no astrophysical applications in strong-field 
general relativistic flows have yet been attempted with this formulation. 



2.2.2 Papadopoulos and Font 

In this formulation [222] the spatial velocity components of the 4- velocity, u*, 
together with the rest-frame density and internal energy, p and £, provide a 
unique description of the state of the fluid at a given time and are taken as 
the primitive variables. They constitute a vector in a five dimensional space 
w = {p,u^,e). The initial value problem for equations (7) and (8) is defined 
in terms of another vector in the same fluid state space, namely the conserved 
variables, U, individually denoted {D,S^,E): 

D = V° = J° = pu°, (40) 

= = r"^ = phv^jj' + pg"' , (41) 
E = U'' = T°° = phu°u°+pg°°. (42) 

Note that the state vector variables slightly differ from previous choices (see, 
e.g., Eqs. (22), (31), (32), (33) and (39)). With those definitions the equations 
of general relativistic hydrodynamics take the standard conservation law form, 

= S, (43) 

with A = (0, z,4). The flux vectors and the source terms S (which depend 
only on the metric, its derivatives and the undifferentiated stress energy tensor), 
are given by 

F^' = ( J^' , T^' , T^°) = {pu^ , phu'u^ + pg'^ , phu\^ + pg°^ ) , (44) 



and 



S = (0,-^/^^^T^^-V^^" T^^). (45) 



The state of the fluid is uniquely described using either vector of variables, 
i.e., either U or w, and each one can be obtained from the other via the defi- 
nitions (40)- (42) and the use of the normalization condition for the 4- velocity, 
g^i.u^u'^ = — 1. The local characteristic structure of the above system of equa- 
tions was presented in [222], where the formulation proved well suited for the 
numerical implementation of HRSC schemes. The formulation presented in this 
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section has been developed for a perfect fluid EOS. Extensions to account for 
generic EOS are given in [220]. This reference further contains a comprehensive 
analysis of general rclativistic hydrodynamics in conservation form. 

A technical remark must be included here: In all conservative formulations 
discussed in Sections 2.1.3, 2.2.1 and 2.2.2, the time update of a given numer- 
ical algorithm is applied to the conserved quantities U. After this update the 
vector of primitive quantities w must be reevaluated, as those are needed in the 
Riemann solver (see Section 3.1.2). The relation between the two sets of vari- 
ables is, in general, not in closed form and, hence, the recovery of the primitive 
variables is done using a root-finding procedure, typically a Newton-Raphson 
algorithm. This feature, distinctive of the equations of (special and) general rel- 
ativistic hydrodynamics - it does not exist in the Newtonian limit - may lead in 
some cases to accuracy losses in regions of low density and small speeds, apart 
from being computationally inefficient. Specific details on this issue for each 
formulation of the equations can be found in Refs. [25, 82, 222]. In particular, 
for the covariant formulation discussed in Section 2.2.1, there exists an analytic 
method to determine the primitive variables, which is, however, computation- 
ally very expensive since it involves many extra variables and solving a quartic 
polynomial. Therefore, iterative methods are still preferred [82]. On the other 
hand, we note that the covariant formulation discussed in this section, when ap- 
plied to null spacctimc foliations, allows for a simple and explicit recovery of the 
primitive variables, as a consequence of the particular form of the Bondi-Sachs 
metric. 

Lightcone hydrodynamics: A comprehensive numerical study of the for- 
mulation of the hydrodynamic equations discussed in this section was presented 
in [222], where it was applied to simulate onc-dimcnsional rclativistic flows on 
null (lightlike) spacetime foliations. The various demonstrations performed in- 
clude standard shock tube tests in Minkowski spacetime, perfect fluid accretion 
onto a Schwarzschild black hole using ingoing null Eddington-Finkclstcin coor- 
dinates, dynamical spacetime evolutions of rclativistic polytropes (i.e., stellar 
models satisfying the so-called Tolman-Oppenheimer-Volkoff equations of hy- 
drostatic equilibrium) sliced along the radial null cones, and accretion of self- 
gravitating matter onto a central black hole. 

Procedures for integrating various forms of the hydrodynamic equations on 
null hypersurfaces are much less common than on spacelike (3+1) hypcrsurfaccs. 
They have been presented before in [136] (see [34] for a recent implementation). 
This approach is geared towards smooth isentropic flows. A Lagrangian method, 
limited to spherical symmetry, was developed by [182]. Recent work in [75] 
includes an Eulerian, non-conservative, formulation for general fluids in null 
hypersurfaces and spherical symmetry, including their matching to a spacelike 
section. 

The general formalism laid out in [222, 220] is currently being systemati- 
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cally applied to astrophysical problems of increasing complexity. Applications 
in spherical symmetry include the investigation of accreting dynamic black holes, 
which can be found in [222, 223]. Studies of the gravitational collapse of super- 
massive stars are discussed in [159] and studies of the interaction of scalar fields 
with relativistic stars are presented in [270] . Axisymmetric neutron star space- 
times have been considered in [269] , as part of a broader program aimed at the 
study of relativistic stellar dynamics and gravitational collapse using characteris- 
tic numerical relativity. We note that there has been already a proof-of-principle 
demonstration of the inclusion of matter fields in three dimensions [34]. 

2.3 Going further: non-ideal hydrodynamics 

Formulations of the equations of non-ideal hydrodynamics in general relativity 
are also available in the literature. The term 'iion-ideal" accounts for additional 
physics such as viscosity, magnetic fields and radiation. These non-adiabatic 
effects can play a major role in some astrophysical systems as, e.g., accretion 
disks or relativistic stars. 

The equations of viscous hydrodynamics, the Navier-Stokes-Fourier equa- 
tions, have been formulated in relativity in terms of causal dissipative relativis- 
tic fluids (see the Living Reviews article by Miiller [193] and references therein). 
These extended fluid theories, however, remain unexplored, numerically, in as- 
trophysical systems. The reason may be the lack of an appropriate formulation 
well-suited for numerical studies. Work in this direction was done by Peitz 
and Appl [225] who provided a 3-1-1 coordinate- free representation of different 
types of dissipative relativistic fluid theories which possess, in principle, the 
potentiality of being well adapted to numerical applications. 

The inclusion of magnetic fields and the development of formulations for 
the MHD equations, attractive to numerical studies, is still very limited in 
general relativity. Numerical approaches in special relativity arc presented 
in [146, 292, 24]. In particular, Komissarov [146], and Balsara [24] developed 
two difii'erent upwind HRSC (or Godunov-type) schemes, providing the charac- 
teristic information of the corresponding system of equations, and proposed a 
battery of tests to validate numerical MHD codes. 3-1-1 representations of rel- 
ativistic MHD can be found in [272, 84]. In [313] the transport of energy and 
angular momentum in magneto-hydrodynamical accretion onto a rotating black 
hole was studied adopting Wilson's formulation for the hydrodynamic equations 
(conveniently modified to account for the magnetic terms), and the magnetic 
induction equation was solved using the constrained transport method of [84]. 
Recently, Koide et al. [144, 145] performed the first MHD simulation, in general 
relativity, of magnetically driven relativistic jets from an accretion disk around 
a Schwarzschild black hole (see Section 4.2.2). These authors used a second- 
order finite difference central scheme with nonlinear dissipation developed by 
Davis [65]. Even though astrophysical applications of Godunov-type schemes 
(see Section 3.1.2) in general relativistic MHD are still absent, it is realistic to 
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believe this situation may change in the near future. 

The interaction between matter and radiation fields, present in different 
levels of complexity in all astrophysical systems, is described by the equations of 
radiation hydrodynamics. The Newtonian framework is highly developed (see, 
e.g. [181]; the special relativistic transfer equation is also considered in that 
reference). Pons et al. [231] discuss a hyperbolic formulation of the radiative 
transfer equations, paying particular attention to the closure relations and to 
extend HRSC schemes to those equations. General relativistic formulations of 
radiative transfer in curved spacetimes are considered in, e.g., [237] and [316] 
(see also references therein). 

3 Numerical schemes 

We turn now to describe the numerical schemes, mainly those based on fi- 
nite differences, specifically designed to solve non-linear hyperbolic systems of 
conservation laws. As discussed in the previous section, the equations of gen- 
eral relativistic hydrodynamics fall in this category, irrespective of the formula- 
tion. Even though we also consider schemes based on artificial viscosity tech- 
niques, the emphasis is given on the so-called high-resolution shock-capturing 
(HRSC) schemes (or Godunov-type methods), based on (either exact or approx- 
imate) solutions of local Riemann problems using the characteristic structure 
of the equations. Such finite difference schemes (or, in general, finite volume 
schemes) have been the subject of diverse review articles and textbooks (see, 
e.g., [155, 156, 288, 131]). For this reason only the most relevant features will 
be covered here, addressing the reader to the appropriate literature for further 
details. In particular, an excellent introduction on the implementation of HRSC 
schemes in special relativistic hydrodynamics is presented in the Living Reviews 
article by Marti and Miiller [169]. Alternative techniques to finite differences, 
such as Smoothed Particle Hydrodynamics, (pseudo-) Spectral Methods and 
others, are briefly considered last. 

3.1 Finite difference schemes 

Any system of equations presented in the previous section can be solved nu- 
merically by replacing the partial derivatives by finite differences on a discrete 
numerical grid and then advancing the solution in time via some time-marching 
algorithm. Hence, specification of the state- vector U on an initial hypersurface, 
together with a suitable choice of EOS, followed by a recovery of the primitive 
variables, leads to the computation of the fluxes and source terms. Through 
this procedure the first time derivative of the data is obtained, which then leads 
to the formal propagation of the solution forward in time, with a time-step 
constrained by the Courant-Priedrichs-Lewy (CFL) condition. 
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The hydrodynamic equations (either in Newtonian physics or in general rel- 
ativity) constitute a non-linear hyperbolic system and, hence, smooth initial 
data can transform into discontinuous data (cross of characteristics in the case 
of shocks) in a finite time during the evolution. As a consequence, classical fi- 
nite difference schemes (see e.g. [155, 288]) present important deficiencies when 
dealing with such systems. Typically, first order accurate schemes are much 
too dissipative across discontinuities (excessive smearing) and second order (or 
higher) schemes produce spurious oscillations near discontinuities which do not 
disappear as the grid is refined. To avoid these effects standard finite difference 
schemes have been conveniently modified in various ways which ensure high- 
order, oscillation-free accurate representations of discontinuous solutions, as we 
discuss next. 

3.1.1 Artificial viscosity approach 

The idea of modifying the hydrodynamic equations by introducing artificial vis- 
cosity terms to damp the amplitude of spurious oscillations near discontinuities 
was originally proposed by von Neumann and Richtmyer [296] in the context of 
the (classical) Euler equations. The basic idea is to introduce a purely artificial 
dissipative mechanism whose form and strength are such that the shock tran- 
sition becomes smooth, extending over a small number of intervals Aa; of the 
space variable. In their original work von Neumann and Richtmyer proposed 
the following expression for the viscosity term: 



with a = p{kAxy^, v being the fluid velocity, p the density. Ax the spatial 
interval, and k a constant parameter whose value is conveniently adjusted in 
every numerical experiment. This parameter controls the number of zones in 
which shock waves are spread. 

This type of generic recipe, with minor modifications, has been used in 
conjuction with standard finite difference schemes in all numerical simulations 
employing May and White's formulation, mostly in the context of gravitational 
collapse, as well as Wilson's formulation. So, for example, in May and White's 
original code [173] the artificial viscosity term, obtained in analogy with the one 
proposed by von Neumann and Richtmyer [296] , is introduced in the equations 
accompanying the pressure, in the form: 



Further examples of similar expressions for the artificial viscosity terms, in 
the context of Wilson formulation, can be found in, e.g., [299, 126]. A state-of- 
the-art formulation of the artificial viscosity approach is reported in [17]. 
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The main advantage of the artificial viscosity approach is its simphcity, which 
results in high computational efficiency. Experience has shown, however, that 
this procedure is both, problem dependent and inaccurate for ultra relativistic 
flows [209, 17]. Furthermore, the artificial viscosity approach has the inherent 
ambiguity of finding the appropriate form for Q that introduces the necessary 
amount of dissipation to reduce the spurious oscillations and, at the same time, 
avoids introducing excessive smearing in the discontinuities. In many instances 
both properties are difficult to achieve simultaneously. A comprehensive nu- 
merical study of artificial viscosity induced errors in strong shock calculations 
in Newtonian hydrodynamics (including also proposed improvements) was pre- 
sented by Noh [208]. 

3.1.2 High-resolution shock-capturing (HRSC) upwind schemes 

In finite difference schemes, convergence properties under grid refinement must 
be enforced to ensure that the numerical results are correct (i.e., if a scheme with 
an order of accuracy a is used, the global error of the numerical solution has 
to tend to zero as 0(Aa;)" as the cell width Aa; tends to zero). For hyperbolic 
systems of conservation laws, schemes written in conservation form are preferred 
since, according to the Lax-Wcndroff theorem [153], they guarantee that the 
convergence, if it exists, is to one of the so-called weak solutions of the original 
system of equations. Such weak solutions are generalized solutions that satisfy 
the integral form of the conservation system. They arc classical solutions 
(continuous and differentiable) in regions where they are continuous and have 
a finite number of discontinuities. 

For the sake of simplicity let us consider in the following an initial value 
problem for a one-dimensional scalar hyperbolic conservation law. 



and let us introduce a discrete numerical grid of space-time points {xj,t"). An 
explicit algorithm written in conservation form updates the solution from time 
to the next time level as: 



p + q + I arguments and At and Ax are the time-step and cell width respec- 
tively. Furthermore, u" is an approximation to the average of u{x, t) within the 
numerical cell [xj_i/2,Xj+i/2] (a;j±i/2 = {xj +Xj±i)/2): 



du ^ df{u) 
dt dx 



u{x,t = 0) = uo{x) 



(46) 




(47) 

u) = f{u)) of 




(48) 
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The class of all weak solutions is too wide in the sense that there is no unique- 
ness for the initial value problem. The numerical method should, hence, guar- 
antee convergence to the physically admissible solution. This is the vanishing- 
viscosity limit solution, that is, the solution when rj ^ 0, of the "viscous version" 
of the scalar conservation law, Eq. (46): 

Mathematically, the solution one needs to search for is characterized by the 
so-called entropy condition (in the language of fluids, the condition that the 
entropy of any fluid element should increase when running into a discontinuity). 
The characterization of these entropy- satisfying solutions for scalar equations 
was given by Olcinik [213]. For hyperbolic systems of conservation laws it was 
developed by Lax [152]. 

The Lax-Wendroff theorem [153] cited above does not establish whether the 
method converges. To guarantee convergence, some form of stability is required, 
as Lax first proposed for linear problems [Lax equivalence theorem; see, e.g., 
[241]). Along this direction, the notion of total-variation stability has proven 
very successful although powerful results have only been obtained for scalar 
conservation laws. The total variation of a solution at time t = t", TV(u"), is 
defined as 

TV(n")=^|t.^+i-t.,"|. (50) 

A numerical scheme is said to be TV-stable if TV(m") is bounded for all At at 
any time for each initial data. In the case of non-linear, scalar conservation laws 
it can be proved that TV-stability is a sufficient condition for convergence [155], 
as long as the numerical schemes are written in conservation form and have 
consistent numerical flux functions. Current research has focused on the devel- 
opment of high-resolution numerical schemes in conservation form satisfying the 
condition of TV-stability, e.g., the so-called total variation diminishing (TVD) 
schemes [119] (see below). 

Let us now consider the specific system of hydrodynamic equations as formu- 
lated in Eq. (35) and let us consider a single computational cell of our discrete 
spacetime. Let be a region (simply connected) of a given four-dimensional 
manifold Ai, bounded by a closed three-dimensional surface dil. We further 
take the 3-surface dfl as the standard-oriented hyper-parallelepiped made up 
of two spacelike surfaces {Yi^o,'E,^o^^^o} plus timelike surfaces {S^,;, Sj-i^^j,;} 
that join the two temporal slices together. By integrating system (35) over a 
domain fl of a given spacetime, the variation in time of the state vector U 
within is given keeping apart the source terms - by the fluxes through 
the boundary dfl. The integral form of system (35) is 

Jn V=5 dx^ Jn 9x' Jq 
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which can be written in the fohowing conservation form, weh-adapted to nu- 
merical appUcations: 

(UAl/),o+AxO - (UAy),o = 

^F^dx°dx^dx^ - [ ^F^dx°dx'^dx^ ] 

-gF^dx^dx^dx"^ - J ^/^F^dx^dx^dx"^^ 
+ I Sdfl, (52) 




where 



U = -3- / y/^Udx^dx'^dx^ 
Jav 



(53) 



Ay 



= / / ^dx^dx'^dx^. (54) 

J J x"^ J x^ 

A numerical scheme written in conservation form ensures that, in the ab- 
sence of sources, the (physically) conserved quantities, according to the partial 
differential equations, are numerically conserved by the finite difference equa- 
tions. 

The computation of the time integrals of the interface fluxes appearing in 
Eq. (52) is one of the distinctive features of upwind HRSC schemes. One needs 
first to define the so-called numerical fluxes, which are recognized as approxi- 
mations to the time-averaged fluxes across the cell-interfaces, which depend on 
the solution at those interfaces, \]{x^ + Ax^ /2,x^), during a time step, 

^^+i'^iiC (55) 

At the cell-interfaces the flow can be discontinuous and, following the sem- 
inal idea of Godmiov [112], the numerical fluxes can be obtained by solving a 
collection of local Riemann problems, as is depicted in Fig. 2. This is the ap- 
proach followed by the so-called Godunov-type methods [121, 79]. Fig. 2 shows 
how the continuous solution is locally averaged on the numerical grid, a process 
which leads to the appearance of discontinuities at the cell-interfaces. Phys- 
ically, every discontinuity decays into three elementary waves: a shock wave, 
a rarefaction wave and a contact discontinuity. The complete structure of the 
Riemann problem can be solved analytically (see [112] for the solution in New- 
tonian hydrodynamics and [167, 232] in special relativistic hydrodynamics) and, 
accordingly, used to update the solution forward in time. 
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Figure 2: Godunov's scheme: local solutions of Riemann problems. At every 

interface, Xj_^_i and x^^s^ a local Riemann problem is set up as a result 

of the discretization process (bottom panel), when approximating the numerical 
solution by piecewise constant data. At time these discontinuities decay into 
three elementary waves which propagate the solution forward to the next time 
level t""*"^ (top panel). The time step of the numerical scheme must satisfy the 
Courant-Friedrichs-Lewy condition, being small enough to prevent the waves 
from advancing more than Ax/ 2 in At. 



For reasons of numerical efficiency and, particularly in multidimensions, the 
exact solution of the Riemann problem is frequently avoided and linearized 
(approximate) Riemann solvers are preferred. These solvers are based on the 
exact solution of Riemann problems corresponding to a linearized version of 
the original system of equations. After extensive experimentation, the results 
achieved with approximate Riemann solvers are comparable to those obtained 
with the exact solver (see [288] for a comprehensive overview of Godunov-type 
methods, and [169] for an excellent summary of approximate Riemann solvers 
in special relativistic hydrodynamics). 

In the frame of the local characteristic approach the numerical fluxes ap- 
pearing in Eq. (52) are computed according to some generic flux-formula which 
makes use of the characteristic information of the system. For example, in Roe's 
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approximate Riemann solver [242] it adopts the following functional form; 

F^-+i = \ I^F(wr) + F(wl) - |An|A5„?„j , (56) 

where wl and wr are the values of the primitive variables at the left and right 
sides, respectively, of a given cell interface. They are obtained from the cell- 
centered quantities after a suitable monotone reconstruction procedure. 

The way these variables are computed determines the spatial order of accu- 
racy of the numerical algorithm and controls the amplitude of the local jumps 
at every cell interface. If these jumps arc nionotonically reduced, the scheme 
provides more accurate initial guesses for the solution of the local Riemann prob- 
lems (the average values give only first-order accuracy). A wide variety of cell 
reconstruction procedures is available in the literature. Among the slope limiter 
procedures (see, e.g., [288, 156]) most commonly used for TVD schemes [119] are 
the second order, piecewise linear reconstruction, introduced by van Leer [291] 
in the design of the MUSCL scheme (Monotonic Upstream Scheme for Conser- 
vation Laws), and the third order, piecewise parabolic reconstruction developed 
by Colella and Woodward [62] in their Piecewise Parabolic Method (PPM). 
Since TVD schemes are only first-order accurate at local extrema, alternative 
reconstruction procedures where some growth of the total variation is allowed 
have also been developed. Among those we mention the total variation hounded 
(TVB) schemes [268] and the essentially nan- oscillatory (ENO) schemes [120]. 

Alternatively, high-order methods for non-linear hyperbolic systems have 
also been designed using flux limiters rather than slope limiters, as in the FCT 
scheme; of Boris and Book [49]. In this approach, the numerical flux consists of 
two pieces, a high-order flux (e.g. the Lax-Wendroff flux) for smooth regions 
of the flow, and a low-order flux (e.g. the flux from some monotone method) 
near discontinuities, F = F;,, - (1 — ^){Ph — Fi) with the limiter $ e [0,1] 
(see [288, 156] for further details). 

The last term in the flux- formula, Eq. (56), represents the numerical viscosity 
of the scheme, and it makes explicit use of the characteristic information of 
the Jacobian matrices of the system. This information is used to provide the 
appropriate amount of numerical dissipation to obtain accurate representations 
of discontinuous solutions without excessive smearing, avoiding, at the same 
time, the growth of spurious numerical oscillations associated with the Gibbs 
phenomenon. In Eq. (56), {A„, r„}„=i...5 are, respectively, the eigenvalues and 
right-eigenvectors of the Jacobian matrix of the flux vector. Correspondingly, 
the quantities {Aa;„}„=i...5 are the jumps of the so-called characteristic variables 
across each characteristic fleld. They are obtained by projecting the jumps of 
the state-vector variables with the left-eigenvectors matrix: 

5 

U(wr) - U(wl) = ^<^nrn. (57) 

n=l 
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The "tilde" in Eqs. (56) and (57) indicates that the corresponding fields are 
averaged at the cell interfaces from the left and right (reconstructed) values. 

During the last few years most of the standard Ricmann solvers developed 
in Newtonian fluid dynamics have been extended to relativistic hydrodynamics; 
Eulderink [82], as discussed in Section 2.2.1, explicitly derived a relativistic 
Roe's Riemann solver [242], Schneider ct al. [250] carried out the extension 
of Harten, Lax, van Leer and Einfeldt's (HLLE) method [121, 79], Martf and 
Miiller [168] extended the PPM method of Woodward and Colella [306], Wen 
et al [297] extended Glimm's exact Riemann solver, Dolezal and Wong [72] put 
into practice Shu-Osher ENO techniques, Balsara [23] extended Colella's two- 
shock approximation and Donat et al. [73] extended Marquina's method [74]. 
Recently, much effort has been spent concerning the exact special relativistic 
Riemann solver and its extension to multidimensions [167, 232, 238, 239]. The 
interested reader is addressed to [169] and references therein for a comprehensive 
description of such solvers in special relativistic hydrodynamics. 

3.1.3 High-order central schemes 

The use of high-order non-oscillatory symmetric (central) TVD schemes for 

solving hyperbolic systems of conservation laws, emerged at the mid 1980s 
[65, 243, 311, 206] (see also [312] and [288] and references therein) as an al- 
ternative approach to upwind HRSC schemes. Symmetric schemes are based 
on either high-order schemes (e.g. Lax-Wendroff) with additional dissipative 
terms [65, 243, 311], or on non-oscillatory high-order extensions of first-order 
central schemes (e.g. Lax-Priedrichs) [206]. One of the nicest properties of cen- 
tral schemes is that they exploit the conservation form of the Lax-Wendroff or 
Lax-Friedrichs schemes. Therefore, they yield the correct propagation speeds 
of all nonlinear waves appearing in the solution. Furthermore, central schemes 
sidestep the use of Ricmann solvers, which results in enhanced computational 
efficiency, especially in multidimensional problems. Its use is, thus, not limited 
to those systems of equations where the characteristic information is explicitly 
known or to systems where the solution of the Ricmann problem is prohibitively 
expensive to compute. This approach has gradually developed during the last 
decade to reach a mature status where a number of straightforward central 
schemes of high order can be applied to any nonlinear hyperbolic system of 
conservation laws. The typical results obtained for the Euler equations show a 
quality comparable to that of upwind HRSC schemes, at the expense of a small 
loss of sharpness of the solution at discontinuities [288] . An up-to-date summary 
of the status and applications of this approach is discussed in [288, 148, 281]. 

In the context of special and general relativistic MHD, Koide et al [144, 145] 
applied a second-order central scheme with nonlinear dissipation developed by 
[65] to the study of black hole accretion and formation of relativistic jets. One- 
dimensional test simulations in special relativistic hydrodynamics performed by 
the author and coworkers [96] using the SLIC (slope limiter centred) scheme (see 
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[288] for details) showed its capabilities to yield satisfactory results, comparable 
to those obtained by HRSC schemes based on Riemann solvers, even well in- 
side the ultrarelativistic regime. The slopes of the original central scheme were 
limited using high-order reconstruction procedures such as PPM [62], which 
was essential to keep the inherent diffusion of central schemes at discontinuities 
at reasonable levels. Very recently, Del Zanna and Bucciantini [67] assessed 
a third-order convex essentially non-oscillatory central scheme in multidimen- 
sional special relativistic hydrodynamics. Again, these authors obtained results 
as accurate as those of upwind HRSC schemes in standard tests (shock tubes, 
shock reflection test). Yet another central scheme has been assessed by [17] in 
one-dimensional special and general relativistic hydrodynamics, where similar 
results to those of [67] are reported. These authors also validate their central 
scheme in simulations of spherical accretion on to a Schwarzschild black hole, 
and further provide a comparison with an artificial viscosity based scheme. 

It is worth to underline that early pioneer approaches in the field of rela- 
tivistic hydrodynamics [209, 58] used standard finite difference schemes with 
artificial viscosity terms to stabilize the solution at discontinuities. However, as 
discussed in Section 2.1.2, those approaches only succeeded in obtaining accu- 
rate results for moderate values of the Lorentz factor {W ~ 2). A key feature 
lacking in those investigations was the use of a conservative approach for both, 
the system of equations and the numerical schemes. A posteriori, and to the 
light of the results reported by [67, 17, 96], it appears that this was the ultimate 
reason preventing the extension of the early computations to the ultrarelativistic 
regime. 

The alternative of using high-order central schemes instead of upwind HRSC 
schemes becomes apparent when the spectral decomposition of the hyperbolic 
system under study is not known. The straightforwardness of a central scheme 
makes its use very appealing, especially in multi-dimensions where computa- 
tional efficiency is an issue. Perhaps the most important example in relativistic 
astrophysics is the system of (general) relativistic MHD equations. Despite some 
recent progress has been accomplished in recent years (see, e.g., [24, 146]), much 
more work is needed concerning their solution with HRSC schemes based upon 
Riemann solvers. Meanwhile, an obvious choice is the use of central schemes 
[144, 145]. 

3.1.4 Source terms 

Most "conservation laws" involve source terms on the right hand side of the 
equations. In hydrodynamics, for instance, those terms arise when considering 
external forces such as gravity, which make the right hand side of the momen- 
tum and energy equations no longer zero (see Section 2). Other effects leading 
to the appearance of source terms are radiative heat transfer (accounted for in 
the energy equation) or ionization (resulting in a collection of non-homogeneous 
continuity equations for the mass of each species, which is not conserved sep- 
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arately). The incorporation of the source terms in the solution procedure is a 
common issue to all numerical schemes considered so far. Since a detailed dis- 
cussion on the numerical treatment of source terms is beyond the scope of this 
article, we simply provide some basic information in this section, addressing the 
interested reader to [288, 156] and references therein for further details. 

There arc, essentially, two basic ways of handling source terms. The first 
approach is based on unsplit methods by which a single finite difference formula 
advances the entire equation over one time step (as in Eq. (52)): 

Ur ^ = U," - ^ (f,+i - F,_ . ) + At (58) 

The temporal order of this basic algorithm can be improved by introducing 
successive sub-steps to perform the time update (e.g. predictor-corrector, Shu 
& Osher's conservative high order Runge-Kutta schemes, etc.) 

Correspondingly, the second approach is based on fractional step or splitting 
methods. The basic idea is to split the equation into different pieces (transport 
-|- sources) and apply appropriate methods for each piece independently. In the 
first-order Godunov splitting, U"+'^ = £^*£j'*U", the operator solves the 
homogeneous PDE in the first step to yield the intermediate value U*. Then, 
in the second step, the operator Cf^ solves the ordinary differential equation 
dfU* = S(U*) to yield the final value U""'"^. In order to achieve second-order 
accuracy (assuming each independent method is second order) a common frac- 
tional step method is the Strang splitting, where U""''"'^ = £f^*''^£^*£f^*^^U". 
Therefore, this method advances half time step the solution for the ODE con- 
taining the source terms, and a full time step the conservation law (the PDE) 
in between each source step. 

We note that in some cases the source terms may become stiff, as in phe- 
nomena which cither occur on much faster time scales than the hydrodynamic 
time-step At, or act over much smaller spatial scales than the grid resolution 
Ax. Stiff source terms may easily lead to numerical difhcultics. The interested 
reader is addressed to [156] (and references therein) for further information on 
various approaches to overcome the problems of stiff source terms. 

3.2 Other techniques 

Two of the most frequently employed alternatives to finite difference schemes 
in numerical hydrodynamics are Smoothed Particle Hydrodynamics (SPH) and, 
to a lesser extent, (pscudo-) Spectral Methods. This section contains a brief de- 
scription of both approaches. In addition, we also mention the Field-Dependent 
Variation method and the Finite Element method. We note, however, that both 
of these approaches have barely been used so far in relativistic hydrodynamics. 
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3.2.1 Smoothed particle hydrodynamics 



The Lagrangian particle method SPH, derived independently by Lucy [161] and 
Gingold and Monaghan [108], has shown successful performance to model fluid 
flows in astrophysics and cosmology. Most studies to date consider Newtonian 
flows and gravity, enhanced with the inclusion of the fluid self-gravity. 

In the SPH method a finite set of extended Lagrangian particles replaces the 
continuum of hydrodynamical variables, the finite extent of the particles being 
determined by a smoothing function (the kernel) containing a characteristic 
length scale h. The main advantage of this method is that it does not require 
a computational grid, avoiding mesh tangling and distortion. Hence, compared 
to grid-based finite volume methods, SPH avoids wasting computational power 
in multidimensional applications, when, e.g., modelling regions containing large 
voids. Experience in Newtonian hydrodynamics shows that SPH produces very 
accurate results with a small number of particles (?« 10*^ or even less). However, 
if more than 10^ particles have to be used for certain problems, and self-gravity 
has to be included, the computational power, which grows as the square of 
the number of particles, may exceed the capabilities of current supercomputers. 
Among the limitations of SPH we note the difficulties in modelling systems with 
extremely different characteristic lengths and the fact that boundary conditions 
usually require a more involved treatment than in finite volume schemes. 

Reviews of the classical SPH equations are abundant in the literature (see, 
e.g., [188, 192] and references therein). The reader is addressed to [192] for a 
summary of comparisons between SPH and HRSC methods. 

Recently, implementations of SPH to handle (special) relativistic (and even 
ultra relativistic) fiows have been developed (see, e.g., [61] and references therein). 
However, SPH has been scarcely applied to simulate relativistic flows in curved 
spacetimes. Relevant references include [140, 149, 150, 271]. 

Following [149] let us describe the implementation of an SPH scheme in 
general relativity. Given a function /(x), its mean smoothed value (/(x)), (x = 
{x,y, z)) can be obtained from 



whore W is the smoothing kernel, h the smoothing length and ^/g' dPx' the vol- 
ume element. The kernel must be differentiable at least once, and the derivative 
should be continuous to avoid large fluctuations in the force felt by a particle. 
Additional considerations for an appropriate election of the smoothing kernel 
can be found in [109] . The kernel is required to satisfy a normalization condition 



which is assured by choosing VF(x,x'; h) — ^(x)ri(w), with ti = |x — x'|//i, ^(x) 
a normalization function and Q.{v) a standard spherical kernel. 




(59) 




(60) 
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The smooth approximation of gradients of scalar functions can be written 

as 



(V/(x))=V(/(x))-(/(x))VlnC(x), 



(61) 



and that corresponding to the divergence of a vector reads 



(V • A(x)) = V • (A(x)) - (A(x)) • InC(x). 



(62) 



Discrete representations of these procedures are obtained after introducing 
the number density distribution of particles (n(x)) = J2a=i ^{^"'^a)/ ^/d-: with 
{xa}o=i,...,Ar the collection of A''-particles where the functions are known. The 
previous representations then read: 



with Q.ab = f^(xa, Xb; /i). These approximations can then be used to derive 
the SPH version of the general relativistic hydrodynamic equations. Explicit 
formulae are reported in [149]. The time evolution of the final system of ODEs 
is performed in [149] using a second-order Runge-Kutta time integrator with 
adaptive time step. As in non-Riemann solver based finite volume schemes, in 
SPH simulations involving the presence of shock waves, artificial viscosity terms 
must be introduced as a viscous pressure term [163]. 

Rcccintly, Sieglcr and Riffert [271] have developed a Lagrangian conservation 
form of the general relativistic hydrodynamic equations for perfect fluids, with 
artificial viscosity, in arbitrary background spacetimes. Within that formulation 
these authors [271] have built a general relativistic SPH code using the stan- 
dard SPH formalism as known from Newtonian fluid dynamics (in contrast to 
previous approaches, e.g., [163, 140, 149]). The conservative character of their 
scheme has allowed the modelling of ultrarelativistic flows including shocks with 
Lorentz factors as large as 1000. 

3.2.2 Spectral methods 

The basic principle underlying spectral methods consists of transforming the 
partial differential equations into a system of ordinary differential equations by 
means of expanding the solution in a series on a complete basis. The mathe- 
matical theory of these schemes is presented in [114, 56]. Spectral methods are 




(63) 




(64) 




(65) 
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particularly well suited to the solution of elliptic and parabolic equations. Good 
results can also be obtained for hyperbolic equations as long as no discontinu- 
ities appear in the solution. When a discontinuity is present some amount of 
artificial viscosity must be added to avoid spurious oscillations. In the specific 
case of relativistic problems, where coupled systems of elliptic equations (i.e., 
the Einstein constraint equations) and hyperbolic equations (i.e., hydrodynam- 
ics) must be solved, an interesting strategy is to use spectral methods to solve 
the elliptic equations and HRSC schemes for the hyperbolic ones. Using such 
combined methods first results have been obtained in one-dimensional super- 
nova collapse simulations, both in the framework of a tensor-scalar theory of 
gravitation [210, 212] and in General Relativity [211]. 

Following [44] we illustrate the main ideas of spectral methods considering 
the quasi-linear onc-dimcnsional scalar equation: 

du d'^u , du „ .„ 

with u = u{t,x) and A a constant parameter. In the linear case (A = 0), and 
assuming the function u to be periodic, spectral methods expand the function 
in a Fourier series: 

oo 

u{x, t) = ^ [dfc (i) cos(27rfca;) + bk {t) sin(27rfca;)] . (67) 

From the numerical point of view, the series is truncated for a suitable value of 
k. Hence, Eq. (66), with A = 0, can be rewritten as 

^ = -fcV(t), ^ = -kHk{t). (68) 

Finding a solution of the original equation is then equivalent to solve an "infi- 
nite" system of ordinary differential equations, where the initial values of the 
coefficients Ofc and bk arc given by the Fourier expansion of u{x,0). 

In the non-linear case, A / 0, spectral methods proceed in a more convoluted 
way: first, the derivative of u is computed in the Fourier space. Then, a step 
back to the configuration space is taken through an inverse Fourier transform. 
Finally, after multiplying du/dx by u in the configuration space the scheme 
returns again to the Fourier space. 

The particular set of trigonometric functions used for the expansion in 
Eq. (67) is chosen because it automatically fulfills the boundary conditions, and 
because a fast transform algorithm is available (the latter is no longer an issue 
for today's computers). If the initial or boundary conditions are not periodic, 
Fourier expansion is no longer useful because of the presence of a Gibbs phe- 
nomenon at the boundaries of the interval. Legendre or Chebyshev polynomials 
are, in this case, the most common base of functions used in the expansions 
(see [114, 56] for a discussion on the different expansions). 
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Extensive numerical applications using (pseudo-) spectral methods have been 
undertaken by the LUTH Relativity Group at the Observatoire de Paris in 
Mcudon. The group is focused in the study of compact objects, as well as 
the associated violent phenomena, gravitational collapse and supernova explo- 
sion. They have developed a fully object-oriented library (based on the C-|— |- 
computcr language) called LORENE [160] to implement (multi domain) spec- 
tral methods within spherical coordinates. A comprehensive summary of ap- 
plications in general relativistic astrophysics is presented in [44]. The most 
recent ones deal with the computation of quasi-equilibrium configurations of 
either synchronized or irrotational binary neutron stars in General Relativity 
[116, 284, 283]. Such initial data are currently being used by fully relativis- 
tic, finite difference time-dependent codes (see Section 3.3.2) to simulate the 
coalescence of binary neutron stars. 



3.2.3 Flow field-dependent variation method 

Richardson and Chung [240] have recently proposed the flow ficld-dcpcndcnt 
variation (FDV) method as a new approach for general relativistic (non ideal) 
hydrodynamics computations. In the FDV method, parameters characteristic 
of the flow field arc computed in order to guide the numerical scheme toward 
a solution. The basic idea is to expand the conservation flow variables into a 
Taylor series in terms of the FDV parameters, which are related to variations of 
physical parameters such as the Lorentz factor, the relativistic Reynolds number 
and the relativistic Froude number. 

The general relativistic hydrodynamic equations are expanded in a special 
form of the Taylor series: 

Un+1 ^ un + ^t^}L_ + ^^l^L_ + o{At% (69) 

= -^+''—8^' 0<«^<1' (71) 

with Sa and S(, denoting the first-order and second-order variation parameters, 
using the above expressions, the time update then reads: 

+ — (^+^.^^j+^(Ai^)- (72) 

Combining the conservation form of the equations and the rearranged Taylor 
series expansion allows to rewrite the time update into standard matrix (resid- 
ual) form, which can then be discretized using either standard finite difference 
or finite element methods [240]. 
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The physical interpretation of the coefficients Sa and Sb is the foundation of 
the FDV method. The first-order parameter Sa is proportional to ai9U"+^/9a;j, 
where a; is the convection Jacobian representing the change of convective mo- 
tion. If the Lorentz factor remains constant in space and time then Sa = 0. 
However, if the Lorentz factor between adjacent zones is large, = 1. Similar 
assessments in terms of the Reynolds number can be provided for the diffusion 
and diffusion gradients, while the Froude number is used for the source term 
Jacobian dS/dTJ. Correspondingly, the second-order FDV parameters are 
chosen to be exponentially proportional to the first-order ones. 

Obviously, the main drawback of the FDV method is the dependence of 
the solution procedure on a large number of problem-dependent parameters, 
a drawback shared to some extent by artificial viscosity schemes. Richardson 
and Chung [240] have implemented the FDV method in a C-|— |- code called 
GRAFSS (General Relativistic Astrophysical Flow and Shock Solver). The test 
problems they report are the special relativistic shock tube (problem 1 in the 
notation of [168]) and the Bondi accretion on to a Schwarzschild black hole. 
While their method yields the correct wave propagation, the numerical solution 
near discontinuities has considerably more diffusion than with upwind HRSC 
schemes. Nevertheless, the generality of the approach is worth considering. 
Applications to non ideal hydrodynamics and relativistic MHD are in progress 
[240]. 

3.2.4 Going further 

The finite element method is the preferred choice to solve multidimensional 

structural engineering problems since the late 1960s [318]. First steps in the 
development of the finite element method for modeling astrophysical flows in 
general relativity are discussed by Meier [177]. The method, designed in a fully 
covariant manner, is valid not only for the hydrodynamic equations but also for 
the Einstein and electromagnetic field equations. The most unique aspect of the 
approach presented in [177] is that the grid can be arbitrarily extended into the 
time dimension. Therefore, the standard finite element method generalizes to 
a four-dimensional covariant technique on a spacetime grid, with the engineer's 
"isoparametric transformation" becoming the generalized Lorentz transform. At 
present, the method has shown its suitability to accurately compute the equi- 
librium stellar structure of Newtonian polytropes, either spherical or rotating. 
The main limitation of the finite element method, as Meier points out [177], 
is that it is generally fully implicit, which results in severe computer time and 
memory limitations for three- and four-dimensional problems. 

3.3 State-of-the-eirt three-dimensional codes 

The most advanced time-dependent, finite-difference, three-dimensional Carte- 
sian codes to solve the system of Einstein and hydrodynamics equations are 
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those developed by Shibata [257] and by the Washington University/NCSA/AEI- 
Golm Numerical Relativity collaboration [100, 93]. The main difference between 
both codes lies in the numerical methods used to solve the relativistic hydrody- 
namic equations, artificial viscosity in Shibata's code [257] and upwind HRSC 
schemes in the code of [100, 93]. We note, however, that very recently Shibata 
has upgraded his code to incorporate HRSC schemes (in particular, a Roe-type 
approximate Riemann solver and piecewise parabolic cell-reconstruction proce- 
dures) [259] . Further 3D codes are currently being developed by a group in the 
U.S. (Duez et al. [77]) and by a E.U. Research Training Network collaboration 

[!]• 



3.3.1 Shibata 

In Shibata's code [257] the hydrodynamics equations arc formulated following 
Wilson's approach (Section 2.1.2) for a conformal-traceless reformulation of the 
spacetime variables (see below). An important difference with respect to the 
original system, Eqs. (23)-(25), is that an equation for the entropy is solved in- 
stead of the energy equation. The hydrodynamic equations are integrated using 
van Leer's [291] second order finite difference scheme with artificial viscosity, 
following the approach of a precursor code developed by [215]. 

The ADM Einstein equations are reformulated into a conformal traceless 
system, an idea originally introduced by Shibata and Nakamura [262] (see 
also [198]) and further developed by Baumgartc and Shapiro [29]. This "BSSN" 
formulation, which shows enhanced long-term stability compared to the original 
ADM system, makes use of a conformal decomposition of the 3-metric, 7ij = 
e~'*'^7ij and the trace-free part of the extrinsic curvature, Aij = Kij — jijK/3, 
with the conformal factor cj) chosen to satisfy e^*^ = 7^/^ = det(7jj)^/^. In this 
formulation, as shown by [29], in addition to the evolution equations for the 
conformal three metric 7,;^ and the; conformal-traceless extrinsic ciirvature vari- 
ables Aij , there are evolution equations for the conformal factor 0, the trace of 
the extrinsic curvature K and the "conformal connection functions" F' . Further 
details can be found in [29, 257]. 

The code uses an approximate maximal slicing condition, which amounts to 
solving a parabolic equation for In a, and a minimal distortion gauge condition 
for the shift vector. It also admits 7r-rotation symmetry around the z-axis, as 
well as plane symmetry with respect to the z = plane, allowing computations 
in a quadrant region. In addition, Shibata has also implemented in the code the 
"cartoon" method proposed by the AEI Numerical Relativity group [10], which 
makes possible axisymmetric computations with a Cartesian grid. "Approxi- 
mate" outgoing boundary conditions are used at the outer boundaries, which 
do not completely eliminate numerical errors due to spurious back reflection of 
gravitational waves^. A staggered Leapfrog method is used for the time evo- 

^We note however that all codes based on the 3-1-1 formalism share this problem since 
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lution of the metric quantities. Correspondingly, the hydrodynamic equations 
are updated using a second-order two-step Runge-Kutta scheme. In each time 

step, the staggered metric quantities needed for the hydrodynamics update are 
properly extrapolated to intermediate time levels to reach the desired order of 
accuracy. 

The code developed by Shibata [257, 259] has been tested in a variety of 
problems, including spherical collapse of dust to a black hole, signalled by the ap- 
pearance of the apparent horizon (comparing ID and 3D simulations), stability 
of spherical stars and computation of the radial oscillation period, quadrupole 
oscillations of perturbed spherical stars and computation of the associated gravi- 
tational radiation, preservation of the rotational profile of (approximate) rapidly 
rotating stars, and the preservation of a co- rotating binary neutron star in a 
quasi-equilibrium state (assuming a conformally flat 3-metric) for more than 
one orbital period. Improvements of the hydro dynamical schemes have been 
considered very recently [259] , in order to tackle problems in which shocks play 
an important role, e.g. stellar core collapse to a neutron star. Shibata's code 
has allowed important breakthroughs in the study of the morphology and dy- 
namics of various general relativistic astrophysical problems, such as black hole 
formation resulting from both, the rotational gravitational collapse of neutron 
stars and supermassive stars, and, perhaps most importantly, the coalescence 
of binary neutron stars, a long-standing problem in numerical relativistic hy- 
drodynamics. These applications are discussed in Section 4. The most recent 
simulations of binary neutron star coalescence [265] have been performed on a 
FACOM VPP5000 computer with typical grid sizes of (505,505,253) in {x, y, z). 

3.3.2 CACTUS/GR_ASTRO 

A three-dimensional general relativistic hydrodynamics code developed by the 
Washington Univcrsity/NCSA/AEI-Golm collaboration for the NASA Neutron 
Star Grand Challenge Project [2] is discussed in Refs. [100, 93]. The code 
is built upon the Cactus Computational Toolkit [55]. A version of this code 
which passed the milestone requirement of the NASA Grand Challenge project, 
was released to the community. This code has been benchmarked at over 140 
GFlop/sec on a 1024 node Cray T3E, with a scaling efficiency of over 95%, 
showing the potential for large scale 3D simulations of realistic astrophysical 
systems. 

The hydrodynamics part of the code uses the conservative formulation dis- 
cussed in Section 2.1.3. A variety of state-of-the-art Riemann solvers are im- 
plemented, including a Roe- like solver [242] and Marquina's flux formula [74]. 
The code uses slope-limiter methods to construct second-order TVD schemes by 

the outer boundaries are located at finite radii. Further work on the development of more 
sophisticated boundary conditions is needed to solve this problem. Alternative solutions are 
to follow the light-cone approach developed by Winicour et al [305] or the conformal formalism 
of Priedrich [103] 
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means of monotonic piecewise linear reconstructions of the cell-centered quanti- 
ties for the computation of the numerical fluxes. The standard "minmod" limiter 
and the "monotonizcd central-difference" limiter [290] are implemented. Both 
schemes provide the desired second-order accuracy for smooth solutions, while 
still satisfying the TVD property. In addition, third-order piecewise parabolic 
(PPM) reconstruction has been recently implemented and used in [279]. 

The Einstein equations are solved using different approaches, (i) the standard 
ADM formalism, (ii) a hyperbolic formulation developed by [43], and (iii) the 
BSSN formulation of [198, 262, 29]. The time evolution of both the ADM and 
the BSSN systems can be performed using several numerical schemes [100, 8, 93]. 
Currently, a Leapfrog (non-staggered in time), and an iterative Crank-Nicholson 
scheme have been coupled to the hydrodynamics solver. The code is designed 
to handle arbitrary shift and lapse conditions, which can be chosen as appropri- 
ate for a given spacetime simulation. The AEI Numerical Relativity group [7] is 
currently focused in developing robust gauge conditions for (vacuum) black hole 
spacetimes (see e.g. [11] and references therein). Hence, all advances accom- 
plished here can be incorporated in future versions of the code for non- vacuum 
spacetime evolutions. Similarly, being a general purpose code, a number of dif- 
ferent boundary conditions can be imposed for either the spacetime variables 
or for the hydrodynamical variables. We refer the reader to [100, 8, 93] for 
additional details. 

The code has been subjected to a series of convergence tests [100], with 
many different combinations of the spacetime and hydrodynamics finite differ- 
encing schemes, demonstrating the consistency of the discrete equations with 
the differential equations. The simulations performed in [100] include, among 
others, the evolution of equilibrium configurations of compact stars (solutions 
to the TOV equations) , and the evolution of relativistically boosted TOV stars 
(w = 0.87c) transversing diagonally across the computational domain, test for 
which an exact solution exists. In [8, 9] the improved stability properties of 
the BSSN formulation of the Einstein equations were compared to the ADM 
system. In particular, in [8] a number of strongly gravitating systems were 
simulated, including weak and strong gravitational waves, black holes, boson 
stars and relativistic stars. While the error growth-rate can be decreased by 
going to higher grid resolutions, the BSSN formulation requires grid resolutions 
higher than the ones needed in the ADM formulation to achieve the same accu- 
racy. Furthermore, it was shown in Ref. [93] that the code successfully passed 
stringent long-term evolution tests, such as the evolution of both, spherical and 
rapidly rotating, stationary stellar configurations, and the formation of an ap- 
parent horizon during the collapse of a relativistic star to a black hole. The 
high accuracy of the hydrodynamical schemes employed has allowed the de- 
tailed investigation of the frequencies of radial, quasi-radial and quadrupolar 
oscillations of relativistic stellar models, and use them as a strong assessment 
of the code. The frequencies obtained have been compared to the frequencies 
computed with perturbative methods and in axisymmetric nonlinear evolutions 
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[92]. In all of the cases considered, the code reproduces these results with ex- 
cellent accuracy and is able to extract the gravitational waveforms that might 
be produced during non-radial stellar pulsations. 

4 Hydrodynamical simulations in relativistic as- 
trophysics 

With the exception of the vacuum two-body problem (i.e. the coalescence of 
two black holes), all realistic astrophysical systems and sources of gravitational 
radiation involve matter. Not surprisingly, the joint integration of the equations 
of motion for matter and geometry was in the minds of theorists from the very 
beginning of numerical relativity. 

Nowadays there is a large body of numerical investigations in the litera- 
ture dealing with hydrodynamical integrations in static background spacetimes. 
Most of those are based on Wilson's Eulerian formulation of the hydrodynamic 
equations and use schemes based on finite differences with some amount of ar- 
tificial viscosity. The use of conservative formulations of the equations, and 
the incorporation of the characteristic information in the design of numerical 
schemes, has started in more recent years. 

On the other hand, time-dependent simulations of self-gravitating flows in 
general relativity, evolving the spacetime dynamically with the Einstein equa- 
tions coupled to a hydrodynamic source, constitute a much smaller sample. 
Although there is much interest in this direction, only the spherically symmet- 
ric case (ID) has been extensively studied. In axisymmetry (2D) fewer attempts 
have been made, with most of them devoted to the study of the gravitational 
collapse of rotating stellar cores, black hole formation and the subsequent emis- 
sion of gravitational radiation. Three-dimensional simulations have only started 
more recently. Much effort is nowadays focused on the study of the coalescence 
of compact neutron star binaries (as well as the vacuum black hole binary coun- 
terpart). These theoretical investigations are driven by the emerging possibility 
of soon detecting gravitational waves with the different experimental efforts 
currently underway. The waveform catalogues resulting from time-dependent 
hydrodynamical simulations may provide some help to data analysis groups, 
since the chances for detection may be enhanced through matched-filtering tech- 
niques. 

In the following we review the status of the numerical investigations in three 
astrophysical scenarios all involving strong gravitational fields and, hence, rela- 
tivistic physics: gravitational collapse, accretion onto black holes, and hydrody- 
namical evolution of neutron stars. Relativistic cosmology, another area where 
fundamental advances have been accomplished through numerical simulations, 
is not considered. The interested reader is addressed to the Living Reviews 
article by Anninos [15] and references therein. 
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4.1 Gravitational collapse 

The study of gravitational collapse of massive stars has been largely pursued 
numerically over the years. This problem was already the main motivation 
when May and White built the first one-dimensional code [173, 174]. Such code 
was conceived to integrate the coupled system of Einstein field equations and 
relativistic hydrodynamics, sidestepping Newtonian approaches. 

By browsing through the literature one realizes that the numerical study 
of gravitational collapse has been neatly split in two main directions since the 
early investigations. First, the computational astrophysics community has tra- 
ditionally focused on the physics of the (type Il/Ib/Ic) Supernova mechanism, 
i.e., collapse of unstable stellar iron cores followed by a bounce and explosion. 
Nowadays, numerical simulations are highly developed, as they include rotation 
and detailed state-of-the-art microphysics (e.g., EOS and sophisticated treat- 
ments for neutrino transport). These studies are currently performed, routinely, 
in multidimensions with advanced HRSC schemes. Progress in this direction 
has been achieved, however, at the expense of necessary approximations in the 
treatment of the hydrodynamics and the gravitational field, by using Newtonian 
physics. In this approach the emission of gravitational radiation is computed 
through the quadrupole formula. For reviews of the current status in this di- 
rection see [192, 138] and references therein. 

On the other hand, the numerical relativity community has been more inter- 
ested, since the beginning, in highlighting those aspects of the collapse problem 
having to do with relativistic gravity, in particular the emission of gravitational 
radiation from non-spherical collapse (see [207] for a recent Living Reviews ar- 
ticle on gravitational radiation from gravitational collapse). Much of the effort 
has focused on developing reliable numerical schemes to solve the gravitational 
field equations and much less emphasis, if any, has been given to the details of 
the microphysics of the collapse. In addition, this approach has mostly consid- 
ered gravitational collapse leading to black hole formation, employing relativis- 
tic hydrodynamics and gravity. Quite surprisingly both approaches have barely 
interacted over the years, except for a handful of simulations in spherical sym- 
metry. Nevertheless, numerical relativity is steadily reaching a state in which 
it is not unrealistic to expect a much closer interaction in the near future (see, 
e.g. [261, 258, 93, 259] and references therein). 

4.1.1 Core collapse supernovae 

At the end of their thermonuclear evolution, massive stars in the (Main Se- 
quence) mass range from 9Mq to 30Mq develop a core composed of iron group 
nuclei which becomes dynamically unstable against gravitational collapse. The 
iron core collapses to a neutron star or a black hole releasing gravitational bind- 
ing energy of the order - 3 x lO'"'-"' erg {M/Mq)'^{R/10 km)-i, which is sufRcient 
to power a supernova explosion. The standard model of type Il/Ib/Ic Super- 
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Figure 3: Schematic profiles of the velocity versus radius at three different 

times during core collapse: at the point of "last good homology" , at bounce and 
at the time when the shock wave has just detached from the inner core. 



novae involves the presence of a strong shock wave formed at the edge between 
the homologous inner core and the outer core, which falls at roughly free-fall 
speed. A schematic illustration of the dynamics of this process is depicted in 
Fig. 3. The shock is produced after the bounce of the inner core when its den- 
sity exceeds nuclear density. Numerical simulations have tried to elucidate if 
this shock is strong enough to propagate outwards through the star's mantle 
and envelope given certain initial conditions for the material in the core (an 
issue subject to important uncertainties in the nuclear EOS), as well as through 
the outer layers of the star. In the accepted scenario which has emerged from 
numerical simulations, the existence of the shock wave together with neutrino 
heating that rc-encrgizes it (in the so-called dclaycd-mcchanism by which the 
stalled prompt supernova shock is reactivated by an increase in the post-shock 
pressure due to neutrino energy deposition [301, 33]), and convective overturn 
which rapidly transports energy into the shocked region [63, 32] (and which 
can lead to large-scale deviations from spherically symmetric explosions), are 
all necessary ingredients for any plausible explosion mechanism (see [138] for a 
review) . 

However, the path to reach such conclusions has been mostly delineated from 
numerical simulations in one spatial dimension. Fully relativistic simulations of 
microphysically detailed core collapse off spherical symmetry arc still absent 
and they might well introduce some modifications. The broad picture described 
above has been demonstrated in multidimensions using sophisticated Newtonian 
models, as is documented in [192]. 
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Spherically symmetric simulations. May and White's formulation and 
their corresponding one-dimensional code formed the basis of most spherically 

symmetric codes built to study the collapse problem. Wilson [298] investigated 
the effect of neutrino transport on stellar collapse concluding that heat con- 
duction by neutrinos does not produce the ejection of material, in contrast to 
previous results [64, 18]. The code solved the coupled set of hydrodynamic 
and Einstein equations, supplemented with a Boltzmann transport equation to 
describe the neutrino flow. Van Riper [294] used a spherically symmetric gen- 
eral rclativistic code to study adiabatic collapse of stellar cores, considering the 
purely hydrodynamical bounce as the preferred explosion mechanism. The im- 
portant role of general rclativistic effects to produce collapses otherwise absent 
in Newtonian simulations was emphasized. Bruenn [52] developed a code in 
which the neutrino transport was handled with an approximation, the multi- 
group flux-limited diffusion. Baron et al. [28] obtained a successful and very 
energetic explosion for a model in which the value of the incompressibility mod- 
ulus of symmetric nuclear matter at zero temperature was particularly small, 
j^sym _ -j^gQ ]y[gY ^j^g precise value, nowadays preferred around 240 MeV [111], 
is still a matter of debate). Mayle et al. [175] computed neutrino spectra from 
stellar collapse to stable neutron stars for different collapse models using, as 
in Ref. [52], multigroup flux-limited diffusion for the neutrino transport. Van 
Riper [295] and Bruenn [53] verified that a softer supranuclear EOS, combined 
with general rclativistic hydrodynamics, increases the chances for a prompt 
explosion. In [249, 248] and [179] the neutrino transport was first handled with- 
out approximation by using a general rclativistic Boltzmann equation. Whereas 
in [249, 248] the neutrino transport is described by moments of the general rela- 
tivistic Boltzmann equation in polar slicing coordinates [26], [179] used maximal 
slicing coordinates. Miralles et al. [185], employing a realistic EOS based upon a 
Hartree-Fock approximation for Skyrme-type nuclear forces for densities above 
nuclear density, found that the explosion was favoured by soft EOS, a result 
qualitatively similar to that of Baron et al. [28], who used a phenomenolog- 
ical EOS. Swesty et al. [280] also focused on the role of the nuclear EOS in 
stellar collapse on prompt timescales. Contrary to previous results they found 
that the dynamics of the shock is almost independent of the nuclear ineom- 
pressibility once the EOS is not unphysically softened as in earlier simulations 
(e.g. [294, 28, 295, 53, 185]). Swesty and coworkers used a finite temperature 
compressible liquid drop model EOS for the nuclear matter component [151]. 
For the shock to propagate promptly to a large radius they found that the EOS 
must be very soft at densities just above nuclear densities, which is inconsistent 
with the 1.44M0 neutron star mass constraint imposed by observations of the 
binary pulsar PSR 1913-1-16. 

Prom the above discussion it is clear that numerical simulations show a 
strong sensitivity of the explosion mechanism to the details of the post-bounce 
evolution: general relativity, the nuclear EOS and the properties of the nascent 
neutron star, the treatment of the neutrino transport, and the neutrino-matter 
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interaction. Recently, state-of-the-art treatments of the neutrino transport have 
been achieved, even in the general relativistic case, by solving the Boltzmann 
equation in connection with the hydrodynamics, instead of using multigroup 
flux- limited diffusion [178, 233, 234, 157]. The results of the few spherically 
symmetric simulations available show, however, that the models do not ex- 
plode, neither in the Newtonian or in the general relativistic case. Therefore, 
computationally expensive multidimensional simulations with Boltzmann neu- 
trino transport, able to account for convective effects, are needed to draw further 
conclusions on the viability of the neutrino-driven explosion mechanism. 

From the numerical point of view, many of the above investigations used 
artificial viscosity techniques to handle shock waves. Together with a detailed 
description of the microphysics, the correct numerical modeling of the shock 
wave is the major issue in stellar collapse. In this context, the use of HRSC 
schemes, specifically designed to capture discontinuities, is much more recent, 
starting in the late 1980s with the Newtonian simulations of Fryxell et al. [107] 
using an Eulerian second order PPM scheme (see [192] for a review of the 
present status). There are only a few relativistic simulations, so far restricted 
to spherical symmetry [132, 244, 212]. 

Marti et al. [165] implemented Glaister's approximate Riemann solver [110] 
in a Lagrangian Newtonian hydrodynamical code. They performed a compar- 
ison of the energetics of a stellar collapse simulation using this HRSC scheme 
and a standard May and White code with artificial viscosity, for the same ini- 
tial model, grid size and EOS. They found that the simulation employing a 
Godunov-type scheme produced larger kinetic energies than that corresponding 
to the artificial viscosity scheme, with a factor of two difference in the maximum 
of the infalling velocity. Motivated by this important difFercncc Janka et al. [139] 
repeated this computation with a different EOS, using a PPM second order 
Godunov-type scheme, disagreeing with Marti et al. [165]. The state-of-the- 
art implementation of the tensorial artificial viscosity terms in [139], together 
with the very fine numerical grids employed (unrealistic for three dimensional 
simulations), could be the reason of the discrepancies. 

As mentioned in Section 2.1.3, Godunov-type methods were first used to 
solve the equations of general relativistic hydrodynamics in [166], where the 
characteristic fields of the one-dimensional (spherically symmetric) system were 
derived. The first astrophysical application was stellar collapse. In [41] the hy- 
drodynamic equations were coupled to the Einstein equations in spherical sym- 
metry. The field equations were formulated as a first-order flux-conservative 
hyperbolic system for a harmonic gauge [42] , somehow "resembling" the hydro- 
dynamic equations. HRSC schemes were applied to both systems simultaneously 
(only coupled through the source terms of the equations). Results for simple 
models of adiabatic collapse can be found in [164, 41, 129]. 

A comprehensive study of adiabatic, one-dimensional, core collapse using 
explicit upwind HRSC schemes was presented in [244] (see [309] for a similar 
computation using implicit schemes). In this investigation the equations for the 
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Figure 4: Stills from a movie showing the animations of a relativistic adi- 
abatic core collapse using HRSC schemes (snapshots of the radial profiles of 
various variables are shown at different times). The simulations are taken from 
Ref. [244]: Velocity (top- left), logarithm of the rest-mass density (top-right), 
gravitational mass (bottom-left), and lapse fimction squared (bottom- right). 
See text for details of the initial model. Visualization by Jose V. Romero. 
(To see the movie, please go to the electronic version of this review article at 
http : //www . livingreviews . org/Articles/Volume3/2000-2f ont). 



hydrodynamics and the geometry are written using radial gauge polar slicing 
(Schwarzschild-type) coordinates. The collapse is modeled with an ideal gas 
EOS with a non-constant adiabatic index, which depends on the density as 
r = Tmin + ??(log p — log pb), where pb is the bounce density and 77 = if p < p6 
and ?7 > otherwise [294]. A set of animations of the simulations presented 
in [244] is included in the four MPEG movies in Fig. 4. They correspond to the 
rather stiff Model B of [244]: Tmin = 1.33, r; = 5 and pb = 2.5 x lO^^ gcm'^. 
The initial model is a white dwarf having a gravitational mass of 1.3862M0. The 
animations show the time evolution of the radial profiles of the following fields: 
velocity (movie 1), logarithm of the rest-mass density (movie 2), gravitational 
mass (movie 3) and the square of the lapse function (movie 4). 

The movies show that the shock is sharply resolved and free of spurious os- 
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Figure 5: Still from a movie showing the animation of the time evolution of 
the entropy in a core collapse supernova explosion [141]. The movie shows the 
evolution within the innermost 3000 km of the star and up to 220 ms after 
core bounce. See text for explanation. Visualization by Konstantinos Kifonidis. 
(To see the movie, please go to the electronic version of this review article at 
http : //www. livingreviews . org/Articles/Volume3/2000-2f ont . ) 



dilations. The radius of the inner core at the time of maximum compression, 
at which the infall velocity is maximum (v = —0.62c), is 6.3 km. The gravita- 
tional mass of the inner core at the time of maximum compression is ~ I.OMq. 
The minimum value the quantity reaches is 0.14 which indicates the highly 
relativistic character of these simulations (at the surface of a typical neutron 
star the value of the lapse function squared is ^ 0.75). 



Axisymmetric Newtonian simulations. Beyond spherical symmetry most 
of the existing simulations of core collapse supernova are Newtonian. Axisym- 
metric investigations have been performed by [191, 40, 45] using a realistic EOS 
and some treatment of weak interaction processes but neglecting neutrino trans- 
port, and by [137, 189, 135, 105, 106] employing some approximate description 
of neutrino transport. In addition, [88, 310, 319] have performed Newtonian 
parameter studies of the axisymmetric collapse of rotating polytropes. 

Among the more detailed multi-dimensional non-relativistic hydrodynamical 
simulations are those performed by the MPA/Garching group (an on-line sample 
can be found at their website [190]). As mentioned before these include advanced 
microphysics and employ accurate HRSC integration schemes. To illustrate the 
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degree of sophistication achieved in Newtonian simulations we present in the 
MPEG movie in Fig. 5 an animation of the early evolution of a core collapse 
supernova explosion up to 220 ms after core bounce and shock formation (only 
an intermediate snapshot at 130 ms is depicted in Fig. 5). The movie shows the 
evolution of the entropy within the innermost 3000 km of the star. 

The initial data used in these calculations is taken from the ISM© pre- 
collapse model of a blue supergiant star of [308] . The computations start 20 ms 
after core bounce from a one-dimensional model of [54] . This model is obtained 
with the one-dimensional general relativistic code mentioned previously [52] 
which includes a detailed treatment of the neutrino physics and a realistic EOS, 
and which is able to follow core collapse, bounce and the associated formation 
of the supernova shock. Because of neutrino cooling and energy losses due to 
the dissociation of iron nuclei the shock initially stalls inside the iron core. 

The movie shows how the stalling shock is "revived" by neutrinos streaming 
from the outer layers of the hot, nascent neutron star in the center. A negative 
entropy gradient builds up between the so called "gain-radius", which is the 
position where cooling of hot gas by neutrino emission switches into net energy 
gain by neutrino absorption, and the shock further out. Convective instabilities 
therefore set in, which are characterized by large rising blobs of neutrino heated 
matter and cool, narrow downflows. The convective energy transport increases 
the efficiency of energy deposition in the post-shock region by transporting 
heated material near the shock and cooler matter near the gain radius where 
the heating is strongest. The freshly heated matter then rises again while the 
shock is distorted by the upward streaming bubbles. The reader is addressed 
to [141] for a more detailed description of a more energetic initial model. 

Axisymmetric relativistic simulations. Previous investigations in general 
relativistic rotational core collapse were mainly concerned with the question of 
black hole formation under idealized conditions (see Section 4.1.2), but none 
of those studies has really addressed the problem of supernova core collapse 
which proceeds from white dwarf densities to neutron star densities, involves 
core bounce, shock formation, and shock propagation. 

Wilson [300] first computed neutron star bounces of F = 2 polytropes, and 
measured the gravitational wave emission. The initial configurations were either 
prolate or slightly aspherical due to numerical errors of an otherwise spherical 
collapse. 

More than twenty years later, and with no other simulations in between, 
the first comprehensive numerical study of relativistic rotational supernova core 
collapse in axisymmctry has been performed by Dimmelmeier et al. [68, 69, 
70, 71], who computed the gravitational radiation emitted in such events. The 
Einstein equations were formulated using the so-called conformally flat metric 
approximation [304]. Correspondingly, the hydrodynamic equations were cast 
as the first-order flux-conservative hyperbolic system described in Section 2.1.3. 
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Details of the methodology and of the numerical code are given in [70] . 

Dimmelmeier et al. [71] have simulated the evolution of 26 models in both 
Newtonian and rclativistic gravity. The initial confignrations arc differentially 
rotating relativistic 4/3-poly tropes in equilibrium, which have a central density 
of 10^*^ g cm~^. Collapse is initiated by decreasing the adiabatic index to some 
prescribed fixed value Fi with 1.28 < Fi < 1.325. The EOS consists of a poly- 
tropic and a thermal part for a more realistic treatment of shock waves. Any 
microphysics like electron captures and neutrino transport is neglected. The 
simulations show that the three different types of rotational supernova core col- 
lapse and gravitational waveforms identified in previous Newtonian simulations 
by Zwerger and Miiller [319] (regular collapse, multiple bounce collapse, and 
rapid collapse) are also present in relativistic gravity. However, rotational core 
collapse with multiple bounces is only possible in a much narrower parameter 
range in general relativity. Relativistic gravity has, furthermore, a qualitative 
impact on the dynamics: If the density increase due to the deeper relativistic 
potential is sufficiently large, a collapse which is stopped by centrifugal forces at 
subnuclear densities (and thus undergoes multiple bounces) in a Newtonian sim- 
ulation, becomes a regular, single bounce collapse in relativistic gravity. Such 
collapse type transitions have important consequences for the maximum gravi- 
tational wave signal amplitude. Moreover, in several of the relativistic models 
discussed in [71], the rotation rate of the compact remnant exceeds the criti- 
cal value where MacLaurin spheroids become secularly, and in some cases even 
dynamically, unstable against triaxial perturbations. 

The MPEG movie contained in Fig. 6 shows the time evolution of a multiple 
bounce model (model A2B4G1 in the notation of Ref. [70]), with a type II 
gravitational wave signal. The left panel shows isocontours of the logarithm 
of the density together with the corresponding velocity field distribution. The 
two panels on the right show the time evolution of the gravitational wave signal 
(top panel) and of the central rest-mass density. In the animation the "camera" 
follows the multiple bounces by zooming in and out. The panels on the right 
show how each burst of gravitational radiation coincides with the time of bounce 
of the collapsing core. The oscillations of the nascent proto-neutron star are 
therefore imprinted on the gravitational waveform. Additional animations of 
the simulations performed by [71] can be found at the MPA's website [190]. 

The relativistic models analyzed by Dimmelmeier ct al. [71] cover almost 
the same range of gravitational wave amplitudes (4 x 10^^^ < h'^^ < 3 x 
10~^° for a source at a distance of 10 kpc) and frequencies (60 Hz < v < 
1000 Hz) as the corresponding Newtonian ones [319]. Averaged over all models, 
the total energy radiated in form of gravitational waves is 8.2 x 10"^ Mqc^ 
in the relativistic case, and 3.6 x 10~* Mqc^ in the Newtonian case. For all 
collapse models which are of the same type in both Newtonian and relativistic 
gravity, the gravitational wave signal is of lower amplitude. If the collapse 
type changes, either weaker or stronger signals are found in the relativistic 
case. Therefore, the prospects for detection of gravitational wave signals from 
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Figure 6: Still from an animation showing the time evolution of a rclativis- 
tic core collapse simulation (model A2B4G1 of [71]). Left: Velocity field and 
isocontours of the density. Right: gravitational waveform (top) and central den- 
sity evolution (bottom). Multiple bounce collapse (fizzler), type II signal. The 
camera follows the multiple bounces. Visualization by Harald Dimmelmeier. 
(To see the movie, please go to the electronic version of this review article at 
http : //www. livingreviews . org/Articles/Volume3/2000-2f ont . ) 



axisymmetric supernova rotational core collapse do not improve when taking 
into account relativistic gravity. The signals are within the sensitivity range 
of the first generation laser interferometer detectors if the source is located 
within the Local Group. An on-line waveform catalogue for all models analyzed 
by Dimmelmeier et al. [71] is available at the MPA's website [190]. Finally, 
we note that a program to extend these simulations to full general relativity 
(relaxing the conformally flat metric assumption) has been very recently started 
by Shibata [259]. 

Three-dimensional simulations. To date, there are no three-dimensional 
relativistic simulations of gravitational collapse, in the context of supernova core 
collapse, yet available. All the existing computations have employed Newtonian 
physics. This situation, however, might change in the near future, as very 
recently the first fully relativistic three-dimensional simulations of gravitational 
collapse leading to black hole formation have been accomplished, for rapidly- 
rotating, supermassive neutron stars [261] and supermassive stars [263] (see 
Section 4.1.2), for the head-on collision of two neutron stars [184], and for the 
coalescence of neutron star binaries [257, 264, 265] (see Section 4.3). 

Concerning Newtonian studies Bonazzola and March [45] performed pioneer 
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three-dimensional simulations, using pseudo-spectral methods, very accurate 
and free of numerical or intrinsic viscosity. They confirmed the low amount 
of energy radiated in gravitational waves regardless of the initial conditions of 
the collapse: axisymmetric, rotating or tidally deformed (see also [191]). This 
result only applies to the pre-bounce phase of the supernova collapse as the 
simulations did not consider shock propagation after bounce. More recently, 
[235, 51] have extended the work of [319] studying, with two different PPM 
hydrodynamics codes, the dynamical growth of non-axisymmetric (bar mode) 
instabilities appearing in rotational post-collapse cores. A relativistic extension 
has been performed by [260] (see Section 4.3.1). 

4.1.2 Black hole formation 

Apart from a few onc-dimcnsional simulations, most numerical studies of general 
relativistic gravitational collapse leading to black hole formation have used Wil- 
son's formulation of the hydrodynamics equations and finite difference schemes 
with artificial viscosity. The use of conservative formulations and HRSC schemes 
to study black hole formation, both in two and three dimensions, has started 
more recently. 

Spherically symmetric simulations. Van Riper [294], using a (Lagrangian) 
May and White code, analyzed the mass division between the formation of a 

neutron star or a black hole after gravitational collapse. For the typical (cold) 
EOS used, the critical state was found to lie between the collapses of 1.95M0 
and 1.96M0 cores. 

In [255] a one-dimensional code based on Wilson's hydrodynamical formu- 
lation was used to simulate a general relativistic (adiabatic) collapse to a black 
hole. The fluid equations were finite differenced in a mixed Eulerian-Lagrangian 
form with the introduction of an arbitrary "grid velocity" to ensure sufficient 
resolution throughout the entire collapse. The Einstein equations were formu- 
lated following the ADM equations. Isotropic coordinates and a maximal time 
slicing condition were used to avoid the physical singularity once the black hole 
forms. Due to the non-dynamical character of the gravitational field in the 
case of spherical symmetry (i.e., the metric variables can be computed at every 
time step solving elliptic equations), the code developed by [255] could follow 
relativistic configurations for many collapse time scales ^ GM/c^ after the 
initial appearance of an event horizon. 

A Lagrangian scheme based on the Misner and Sharp [186] formulation for 
spherically symmetric gravitational collapse (as in Ref. [294]) was developed 
by Miller and Motta [182] and later by Baumgarte et al. [30]. The novelty of 
these codes was the use of an outgoing null coordinate u (an "observer-time" 
coordinate, as suggested previously by [127]) instead of the usual Schwarzshild 
time t appearing in the Misner and Sharp equations. Outgoing null coordi- 
nates are ideally suited to study black hole formation as they never cross the 
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black hole event horizon. In these codes the Hernandez-Misner equations [127] 
(or, alternatively, the Misner-Sharp equations) were solved by an explicit finite 
difference scheme similar to the one used by [294]. In [182] the collapse of an 
unstable polytrope to a black hole was first achieved using null slicing. In [30] 
the collapse of a 1.4M0 polytrope with F = 4/3 was compared to the result of 
[294], using the Misner-Sharp equations, finding a 10% agreement. This work 
showed numerically that the use of a retarded time coordinate allows for stable 
evolutions after the black hole has formed. The Lagrangian character of both 
codes has prevented their multidimensional extension. 

Linke et al. [159] analyzed the gravitational collapse of supermassive stars 
in the range 5 x lO^M© - lO^M©. In the same spirit as in Refs. [182, 30], the 
coupled system of Einstein and fluid equations was solved employing coordinates 
adapted to a foliation of the spacetime with outgoing null hypersurfaces. From 
the computed neutrino luminosities, estimates of the energy deposition by uD- 
annihilation were obtained. Only a small fraction of this energy is deposited near 
the surface of the star, where it could cause the ultrarelativistic flow believed 
to be responsible for GRBs. However, the simulations show that for collapsing 
supermassive stars with masses larger than 5 x IQ^Mq, the energy deposition is 
at least two orders of magnitude too small to explain the energetics of observed 
long-duration bursts at cosmological redshifts. 

The interaction of massless scalar fields with sclf-gravitating relativistic 
stars was analyzed in [270] by means of fully dynamic numerical simulations 
of the Einstein-Klein-Gordon perfect fluid system. A sequence of stable, self- 
gravitating, K = 100, iV = 1 relativistic polytropes, increasing the central 
density from pc = 1.5 x 10"^ to 3.0 x 10~^ (G = c = M© = 1) was built. Using 
a compactified spacetime foliation with outgoing null cones, Siebel et al. [270] 
studied the fate of the relativistic stars when they arc hit by a strong scalar 
field pulse with a mass of the order of 10^"^ M©, finding that the star is either 
forced to oscillate in its radial modes of pulsation or to collapse to a black hole 
on a dynamical timcscale. The radiative signals, read off at future null infinity, 
consist of several quasi-normal oscillations and a late time power-law tail, in 
agreement with the results predicted by (linear) perturbation analysis of wave 
propagation in an exterior Schwarzschild geometry. 

Axisymmetric simulations. Beyond spherical symmetry the investigations 

of black hole formation started with the work of Nakamura [194], who first sim- 
ulated general relativistic rotating stellar collapse. He adopted the (2-|-l)-|-l 
formulation of the Einstein equations in cylindrical coordinates and introduced 
regularity conditions to avoid divergence behavior at coordinate singularities 
(the plane z = 0) [196]. The equations were finite differenced using the donor 
cell scheme plus Friedrichs-Lax type viscosity terms. Nakamura used a "hy- 
pergeometric" slicing condition (which reduces to maximal slicing in spherical 
symmetry), which prevents the grid points from hitting the singularity if a black 
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hole forms. The simulations could track the evolution of the collapse of a lOM© 
"core" of a massive star with different amounts of rotational energy and an ini- 
tial central density of 3 x lO""^"^ g cm~^, up to the formation of a rotating black 
hole. However, the resolution affordable due to limitations in computer resources 
(42 X 42 grid points) was not high enough to compute the emitted gravitational 
radiation. Note that the energy emitted in gravitational waves is very small 
compared to the total rest mass energy, which makes its accurate numerical 
computation very challenging. In subsequent works, Nakamura [195] (see also 
[198]) considered a configuration consisting of a neutron star (M = 1.09Mq, 
Pc = 10^^ g cm~'^) with an accreted envelope of O.8IM0, which was thought 
to mimic mass fall-back in a supernova explosion. Rotation and infall velocity 
were added to such configuration, simulating the evolution depending on the 
prescribed rotation rates and rotation laws. 

Bardeen, Stark and Piran, in a series of papers [26, 276, 229, 277], stud- 
ied the collapse of rotating relativistic (F = 2) polytropic stars to black holes, 
succeeding in computing the associated gravitational radiation. The field and 
hydrodynamic equations were formulated using the 3-1-1 formalism, with ra- 
dial gauge and a mixture of (singularity avoiding) polar and maximal slicing. 
The initial model was a spherically symmetric relativistic polytrope in equilib- 
rium of mass M, central density 1.9 x 10^^ {M / Mq)''^ , and radius 6GM/c^ = 
8.8 X lO^M/M© cm. Rotational collapse was induced by lowering the pressure 
in the initial model by a prescribed fraction, and by simultaneously adding an 
angular momentum distribution approximating rigid-body rotation. Their pa- 
rameter space siirvey showed how black hole formation (of the Kerr class) occiirs 
only for angular momenta less than a critical value. The numerical resiilts also 
demonstrated that the gravitational wave emission from axisymmetric rotating 
collapse to a black hole was AE/Mc^ < 7x 10"'*, and that the general waveform 
shape was nearly independent of the details of the collapse. 

Evans [83], building on previous work by Dykema [78], also studied the 
axisymmetric gravitational collapse problem for non-rotating matter configura- 
tions. His numerical scheme to integrate the matter fields was more sophisti- 
cated than in previous approaches, as it included monotonic upwind reconstruc- 
tion procedures and flux limiters. Discontinuities appearing in the flow solution 
were stabilized by adding artificial viscosity terms in the equations, following 
Wilson's approach. 

Most of the axisymmetric codes discussed so far adopted spherical polar co- 
ordinates. Numerical instabilities are encountered at the origin (r = 0) and 
at the polar axis {6 = 0, tt) where some fields diverge due to the coordinate 
singularities. Evans did important contributions towards regularizing the gravi- 
tational field equations in such situations [83] . These coordinate problems have 
deterred researchers from building three-dimensional codes in spherical coordi- 
nates. In this case Cartesian coordinates are adopted which, despite the desired 
property of being everywhere regular, present the important drawback of not 
being adapted to the topology of astrophysical systems. As a result this has 



48 



important consequences on the grid resolution requirements. The only exten- 
sion of an axisymmetric 2+1 code in spherical coordinates to three dimensions 
has been accomplished by Stark [275], although no applications in relativistic 
astrophysics have yet been presented. 

Recently, Alcubierre et al. [10] proposed a method ("cartoon") which does 
not suffer from stability problems at coordinate singularities and where, in 
essence, Cartesian coordinates are used even for axisymmetric systems. Using 
this method, Shibata [258] investigated the effects of rotation on the criterion 
for prompt adiabatic collapse of rigidly and differentially rotating (F = 2) poly- 
tropes to a black hole. Collapse of the initial approximate equilibrium models 
(computed by assuming a conformally flat spatial metric) was induced by a pres- 
sure reduction. In [258] it was shown that the criterion for black hole formation 
depends strongly on the amount of angular momentum, but only weakly on its 
(initial) distribution. Shibata also studied the effects of shock heating using a 
gamma-law EOS, concluding that it is important in preventing prompt collapse 
to black holes in the case of large rotation rates. Using the same numerical 
approach, Shibata and Shapiro[263] have recently studied the collapse of a uni- 
formly rotating supermassive star in general relativity. The simulations show 
that the star, initially rotating at the mass-shedding limit, collapses to form a 
supermassive Kerr black hole with a spin parameter of ~ 0.75. Roughly 90% of 
the mass of the system is contained in the final black hole, while the remaining 
matter forms a disk orbiting around the hole. 

Alternatively, existing ax;isymmetric codes employing the characteristic for- 
mulation of the Einstein equations [305], such as the (vacuum) code presented 
in [113], do not share the axis instability problems of most 2-|-l codes. Hence, 
axisymmetric characteristic codes, once conveniently coupled to hydrodynamics 
codes, are a promising way of studying the axisymmetric collapse problem. First 
steps in this direction are reported in [269], where an axisymmetric Einstein- 
perfect fluid code is presented. This code achieves global energy conservation 
for a strongly perturbed neutron star spacctimc, for which the total energy ra- 
diated away by gravitational waves corresponds to a significant fraction of the 
Bondi mass. 

Three-dimensional simulations. Hydrodynamical simulations of the col- 
lapse of supermassive (F = 2) uniformly rotating neutron stars to rotating black 

holes, using the code discussed in Section 3.3.1, arc presented in [261]. The sim- 
ulations show no evidence for massive disk formation or outflows, which can be 
related to the moderate initial compactness of the stellar models {R/M ~ 6). 
A proof-of-principle of the capabilities of the code discussed in Section 3.3.2 
to study black hole formation is presented in [93], where the gravitational col- 
lapse of a spherical unstable relativistic polytrope is discussed. Similar tests 
with differentially rotating polytropes are given in [77] for a recent artificial 
viscosity-based, three-dimensional general relativistic hydrodynamics code. 
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4.1.3 Critical collapse 

Critical phenomena in gravitational collapse were first discovered numerically 
by Choptuik in spherically-symmetric simulations of the massless Klein-Gordon 
field minimally coupled to gravity [60] . Since then critical phenomena arising at 
the threshold of black hole formation have been found in a variety of physical 
systems, including the perfect fluid model (see [118] for a review). 

Evans and Coleman [85] first observed critical phenomena in spherical col- 
lapse of radiation fluid, i.e. a fluid obeying an EOS p = {T — l)p(l + e) with 
r = 4/3 and £ » 1. The threshold of black hole formation was found for a 
critical exponent T ^ 0.36, in close agreement with that obtained in scalar field 
collapse [60]. Their study used Schwarzschild coordinates in radial gauge and 
polar slicing, and the hydrodynamic equations followed Wilson's formulation. 
Subsequently, Maison [162], using a linear stability analysis of the critical so- 
lution, showed that the critical exponent varies strongly with the parameter 
K = r — 1 of the EOS. More recently, Neilsen and Choptuik [204], using a 
conservative form of the hydrodynamic equations and HRSC schemes, revisited 
this problem. The use of a conservative formulation and numerical schemes 
well adapted to describe ultrarelativistic flows [205], made it possible to find 
evidence of the existence of critical solutions for large values of the adiabatic 
index F (1.89 < F < 2), extending the previous upper limit. 

4.2 Accretion on to black holes 

The study of relativistic accretion and black hole astrophysics is a very ac- 
tive field of research, both theoretically and observationally (see, e.g. [35] and 
references there in). On the one hand, advances in satellite instrumentation, 
e.g., the Rossi X-Ray Timing Explorer (RXTE), and the Advanced Satellite 
for Cosmology and Astrophysics (ASCA), are greatly stimulating and guiding 
theoretical research on accretion physics. The discovery of kHz quasi-periodic 
oscillations in low-mass X-ray binaries, extends the frequency range over which 
these oscillations occur into timescales associated with the innermost regions 
of the accretion process (for a review see [289]). Moreover, in extragalactic 
sources, spectroscopic evidence (broad iron emission lines) increasingly points 
to (rotating) black holes being the accreting central objects [282, 147, 50]. Thick 
accretion discs (or tori) are probably orbiting the central black holes of many 
astrophysical objects such as quasars and other active galactic nuclei (AGNs), 
some X-ray binaries, and microquasars. In addition, they are believed to exist 
at the central engine of cosmic GRBs. 

Disk accretion theory is primarily based on the study of (viscous) station- 
ary flows and their stability properties through linearized perturbations thereof. 
A well-known example is the solution consisting of isentropic constant angular 
momentum tori, unstable to a variety of non-axisymmetric global modes, dis- 
covered by Papaloizou and Pringle [224] (see [22] for a review of instabilities 
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in astrophysical accretion disks). Since the pioneering work by Shakura and 
Sunyaev [253], thin disk models, parametrized by the so-caUed a— viscosity, in 
which the gas rotates with Keplcrian angular momentum which is transported 
radially by viscous stress, have been applied successfully to many astronomical 
objects. The thin disk model, however, is not valid for high luminosity systems, 
as it is unstable at high mass accretion rates. In this regime Abramowicz et 
al. [6] found the so-called slim disk solution, which is stable against viscous and 
thermal instabilities. More recently, much theoretical work has been devoted 
to the problem of slow accretion, motivated by the discovery that many galac- 
tic nuclei are under-luminous (e.g., NGC 4258). Proposed accretion models 
involve the existence of advection-dominated accretion flows (ADAF solution; 
sec, e.g. [203, 201]) and advection-dominated inflow outflow solutions (ADIOS 
solution [36]). The importance of convection for low values of the viscosity pa- 
rameter a is currently being discussed in the so-called convection-dominated 
accretion flow (CDAF solution; see [133] and references therein). The impor- 
tance of magnetic fields and its consequences on the stability properties of this 
solution is critically discussed in [21]. 

For a wide range of accretion problems, the Newtonian theory of gravity 
is adequate for the description of the background gravitational forces (see, 
e.g., [102]). The extensive experience with Newtonian astrophysics has shown 
that explorations of the rclativistic regime benefit from the use of model po- 
tentials. In particular, the Paczynski-Wiita pseudo-Newtonian potential for a 
Schwarzschild black hole [218], gives approximations of general relativistic ef- 
fects with accuracy of 10% — 20% outside the marginally stable radius, which 
corresponds to three times the Schwarzschild radius. Nevertheless, for com- 
prehensive numerical work, a three-dimensional formalism is required, able to 
cover also the maximally rotating hole. In rotating spacetimcs the gravitational 
forces cannot be captured fully with scalar potential formalisms. Additionally, 
geometric regions such as the ergo-spherc would be very hard to model with- 
out a metric description. Whereas the bulk of emission occurs in regions with 
almost Newtonian fields, only the observable features attributed to the inner 
region may crucially depend on the nature of the spacetime. 

In the following we present a summary of representative time- dependent 
accretion simulations in relativistic hydrodynamics. We concentrate on multi- 
dimensional simulations. In the limit of spherical accretion, exact stationary 
solutions exist for both Newtonian gravity [46] and relativistic gravity [ISO]. 
Such solutions are commonly used to calibrate time-dependent hydrodynami- 
cal codes, by analyzing whether stationarity is maintained during a numerical 
evolution [126, 166, 82, 244, 25]. 

4.2.1 Disk accretion 

Pioneering numerical efl[orts in the study of black hole accretion [299, 126, 123, 
124] made use of the the so-called frozen star paradigm of a black hole. In 
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this framework, the time "shcing" of the spacetime is synchronized with that of 
asymptotic observers far from the hole. Within this approach Wilson [299] first 
investigated numerically the timc-dcpcndcnt accretion of inviscid matter onto 
a rotating black hole. This was the first problem to which his formulation of 
the hydrodynamic equations, as presented in Section 2.1.2, was applied. Wilson 
used an axisymmctric hydrodynamical code in cylindrical coordinates to study 
the formation of shock waves and the X-ray emission in the strong-field regions 
close to the black hole horizon, being able to follow the formation of thick 
accretion disks during the simulations. 

Wilson's formulation has been extensively used in time-dependent numer- 
ical simulations of thick disk accretion. In a system formed by a black hole 
surrounded by a thick disk, the gas flows in an effective (gravitational plus cen- 
trifugal) potential, whose structure is comparable with that of a close binary. 
The Roche torus surrounding the black hole has a cusp-like inner edge located 
at the Lagrange point Li where mass transfer driven by the radial pressure gra- 
dient is possible [4]. In [126] (see also [123]) Hawley and collaborators studied, 
in the test-fluid approximation and axisymmetry, the evolution and develop- 
ment of non-linear instabilities in pressure-supported accretion disks formed as 
a consequence of the spiraling infall of fluid with some amount of angular mo- 
mentum. The code used explicit second-order flnite difference schemes with a 
variety of choices to integrate the transport terms of the equations (i.e. those in- 
volving changes in the state of the fluid at a speciflc volume in space) . The code 
also used a staggered grid (with scalars located at the cell centers and vectors 
at the cell boundaries) for its suitability to difference the transport equations. 
Discontinuous solutions were stabilized with artificial viscosity terms. 

With a three-dimensional extension of the axisymmctric code of Hawley et 
al. [125, 126] Hawley [124] studied the global hydrodynamic non-axisymmetric 
instabilities in thick, constant angular momentum accretion gas tori orbiting 
around a Schwarzschild black hole. Such simulations showed that the Papaloizu- 
Pringle instability saturates in a strong spiral pressure wave, not in turbulence. 
In addition, the simulations confirmed that accretion flows through the torus 
could reduce and even halt the growth of the global instability. Extensions to 
Kerr spacctimcs have been recently reported in [66]. 

Yokosawa [313, 314], also using Wilson's formulation, studied the struc- 
ture and dynamics of relativistic accretion disks and the transport of energy 
and angular momentum in magneto-hydrodynamical accretion on to a rotating 
black hole. In his code the hydrodynamic equations are solved using the Flux- 
Corrected- Transport (FCT) scheme [49] (a second-order flux-limiter method 
which avoids oscillations near discontinuities by reducing the magnitude of the 
numerical flux), and the magnetic induction equation is solved using the con- 
strained transport method [84]. The code contains a parametrized viscosity 
based on the a-model [253]. The simulations revealed different flow patterns, 
inside the marginally stable orbit, depending on the thickness, h, of the accre- 
tion disk. For thick disks with h > 4r/j, rn being the radius of the event horizon, 
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Figure 7: Runaway instability of an unstable thick disk: contour levels of the 
rest-mass density p plotted at irregular times from t = 0tot = 11.80 tovh once 
the disk has almost been entirely destroyed. See [91] for details. 

the flow pattern becomes turbulent. 

Igumenshchev and Beloborodov [134] have performed two-dimensional rel- 
ativistic hydrodynamical simulations of inviscid transonic disk accretion on to 
a Kerr black hole. The hydrodynamical equations follow Wilson's formulation 
but the code avoids the use of artificial viscosity. The advection terms are eval- 
uated with an upwind algorithm which incorporates the PPM scheme [62] to 
compute the fluxes. Their numerical work confirms analytic expectations: (i) 
the structure of the innermost disk region strongly depends on the black hole 
spin, and (ii) the mass accretion rate goes as M cx {/\W)^f'^'^^^\ AT-F being 
the potential barrier between the inner edge of the disk and the cusp, and F the 
adiabatic index. 

Thick accretion disks orbiting black holes, on the other hand, may be sub- 
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jected to the so-called runaway instability, as first proposed by Abramowicz et 
al. [5]. Starting from an initial disk filling its Roche lobe, so that mass transfer 
is possible through the cusp located at the Li Lagrange point, two evolutions 
are feasible when the mass of the black hole increases: (i) either the cusp moves 
inwards towards the black hole, which slows down the mass transfer, resulting 
in a stable situation, or (ii) the cusp moves deeper inside the disk material. In 
this case the mass transfer speeds up, leading to the runaway instability. This 
instability, whose existence is still a matter of debate (see e.g. [91] and references 
therein), is an important issue for most models of cosmic GRBs, where the cen- 
tral engine responsible for the initial energy release is such a system consisting 
of a thick disk surrounding a black hole. 

In [91] Font and Daigne presented time-dependent simulations of the run- 
away instability of constant angular momentum thick disks around black holes. 
The study was performed using a fully relativistic hydrodynamics code based on 
HRSC schemes and the conservative formulation discussed in Section 2.1.3. The 
self-gravity of the disk was neglected and the evolution of the central black hole 
was assumed to be that of a sequence of Schwarzschild black holes of varying 
mass. The black hole mass increase is determined by the mass accretion rate 
across the event horizon. In agreement with previous studies based on station- 
ary models, [91] found that by allowing the mass of the black hole to grow the 
disk becomes unstable. For all disk-to-hole mass ratios considered (between 1 
and 0.05), the runaway instability appears very fast on a dynamical timescale 
of a few orbital periods (1 — » 100), typically a few 10 ms and never exceeding 1 
s, for a 2.5 Mq black hole and a large range of mass fluxes (rh > lO^"* Mo/s). 

An example of the simulations performed by [91] appears in Fig. 7. This 
figure shows eight snapshots of the time-evolution of the rest-mass density, from 
t = 0tot = 11.8 torb- The contour levels are linearly spaced with Ap = 0.1p°, 
where p° is the maximum value of the density at the center of the initial disk. 
In Fig. 7 one can clearly follow the transition from a quasi-stationary accretion 
regime (panels (1) to (5)) to the rapid development of the nmaway instability 
in about two orbital periods (panels (6) to (8)). At t = 11.80 torh, the disk 
has almost entirely disappeared inside the black hole whose size has, in turn, 
noticeably grown. 

Extensions of this work to marginally stable (or even stable) constant an- 
gular momentum disks are reported in Zanotti et al. [239] (animations can be 
visualized at the web site listed in Ref. [287]). Furthermore, recent simulations 
with non-constant angular momentum disks and rotating black holes [90] , show 
that the instability is strongly suppressed when adding a small increase out- 
wards of the specific angular momentum of the disk (much smaller than the 
Keplerian value). 



54 



4.2.2 Jet formation 



Numerical simulations of relativistic jets propagating through progenitor stellar 
models of collapsars have been presented in [13]. The coUapsar scenario, pro- 
posed by [307] , is currently the most favoured model for explaining long duration 
GRBs. The simulations performed by [13] employ the three-dimensional code 
GENESIS [12], with a 2D spherical grid and equatorial plane symmetry. The 
gravitational field of the black hole is described by the Schwarzschild metric, and 
the relativistic hydrodynamic equations are solved in the test fluid approxima- 
tion using a Godunov-type scheme. Aloy et al. [13] showed that the jet, initially 
formed by an ad hoc energy deposition of a few 10"'''^ erg s~^ within a 30° cone 
around the rotation axis, reaches the surface of the coUapsar progenitor intact, 
with a maximum Lorentz factor of ^ 33. 

The most promising processes for producing relativistic jets as those observed 
in AGNs, microquasars and GRBs, involve the hydromagnetic centrifugal accel- 
eration of material from the accretion disk [37], or the extraction of rotational 
energy from the ergosphere of a Kerr black hole [226, 39]. Koide and coworkers 
have performed the first MHD simulations of jet formation in general relativity 
[144, 143, 145]. Their code uses the 3+1 formalism of general relativistic conser- 
vation laws of particle number, momentum and energy, and Maxwell equations 
with infinite electric conductivity. The MHD equations are numerically solved 
in the test-fluid approximation (in the background geometry of Kerr spacetime) 
using a finite difference symmetric scheme [65]. The Kerr metric is described in 
Boyer-Lindquist coordinates, with a radial tortoise coordinate to enhance the 
resolution near the horizon. 

In [145], the general relativistic magneto- hydrodynamic behaviour of a plasma 
flowing into a rapidly rotating black hole (a = 0.99995) in a large-scale mag- 
netic fleld is investigated numerically. The initial magnetic fleld is uniform and 
strong, Bq = lOpoc^, Po being the mass density. The initial Alfven speed is 
va = 0.953c. The simulation shows how the rotating black hole drags the iner- 
tial frames around in the ergosphere. The azimuthal component of the magnetic 
fleld increases because of the azimuthal twisting of the magnetic field lines, as is 
depicted in Fig. 8. This frame-dragging dynamo amplifies the magnetic field to 
a value which, by the end of the simulation, is three times larger than the initial 
one. The twist of the magnetic field lines propagates outwards as a torsional 
Alfven wave. The magnetic tension torques the plasma inside the ergosphere 
in a direction opposite to that of the black hole rotation. Therefore, the an- 
gular momentum of the plasma outside receives a net increase. Even though 
the plasma falls into the black hole, electromagnetic energy is ejected along the 
magnetic field lines from the ergosphere, due to the propagation of the Alfven 
wave. By total energy conservation arguments, Koide et al. [145] conclude that 
the ultimate result of the generation of an outward Alfven wave is the magnetic 
extraction of rotational energy of the Kerr black hole, a process the authors call 
MHD Penrose process. Koide and coworkers argue that such process can be 
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Figure 8: Jet formation: twisting of magnetic field lines around a Kerr black 
hole (black sphere) The yellow surface is the ergosphere. The red tubes show 
the magnetic field lines that cross into the ergosphere. Figure taken from [f45] 
(used with permission). 

applicable to jet formation, both in AGNs and GRBs. 

We note that, recently, van Putten and Levinson [293] have considered, 
theoretically, the evolution of an accretion torus in suspended accretion, in con- 
nection with GRBs. These authors claim that the formation of baryon-poor 
outflows may be associated with a dissipative gap in a differentially rotating 
magnetic flux tube supported by an equilibrium magnetic moment of the black 
hole. Numerical simulations of non-ideal MHD, incorporating radiative pro- 
cesses, are necessary to validate this picture. 

4.2.3 Wind accretion 

The term "wind" or hydrodynamic accretion refers to the capture of matter 
by a moving object under the effects of the underlying gravitational field. The 
canonical astrophysical scenario in which matter is accreted in such a non- 
spherical way was suggested originally by Bondi, Hoyle and Lyttleton [128, 47], 
who studied, using Newtonian gravity, the accretion on to a gravitating point 
mass moving with constant velocity through a non-relativistic gas of uniform 
density. The matter flow inside the accretion radius, after being decelerated by 
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a conical shock, is ultimately captured by the central object. Such a process 
applies to describe mass transfer and accretion in compact X-ray binaries, in 
particular in the case in which the donor (giant) star lies inside its Roche lobe 
and loses mass via a stellar wind. This wind impacts on the orbiting compact 
star forming a bow-shaped shock front around it. This process is also believed 
to occur during the common envelope phase in the evolution of a binary system. 

Numerical simulations have extended the simplified analytic models and have 
helped to develop a thorough understanding of the hydrodynamic accretion sce- 
nario, in its fully three-dimensional character (see, e.g., [245, 31] and references 
therein). The numerical investigations have revealed the formation of accretion 
disks and the appearance of non-trivial phenomena such as shock waves and 
tangential {flip-flop) instabilities. 

Most of the existing numerical work has used Newtonian hydrodynamics 
to study the accretion onto non-relativistic stars [245]. For compact accretors 
such as neutron stars or black holes, the correct numerical modeling requires a 
general relativistic hydrodynamical description. Within the relativistic, frozen 
star framework, wind accretion onto "moving" black holes was first studied in 
axisymmetry by Petrich et al. [227]. In this work Wilson's formulation of the 
hydrodynamic equations was adopted. The integration algorithm was borrowed 
from [277] with the transport terms finite-differenced following the prescription 
given in [126]. An artificial viscosity term of the form Q = ap(Av)'^, with a a 
constant, was added to the pressure terms of the equations in order to stabilize 
the numerical scheme in regions of sharp pressure gradients. 

An extensive survey of the morphology and dynamics of relativistic wind 
accretion past a Schwarzschild black hole was later performed by [95, 94]. This 
investigation differs from [227] in both, the use of a conservative formulation for 
the hydrodynamic equations (sec Section 2.1.3) and the use of advanced HRSC 
schemes. Axisymmctric computations were compared to [227] finding major 
differences in the shock location, opening angle and accretion rates of mass and 
momentum. The reasons for the discrepancies are related to the use of different 
formulations, numerical schemes and grid resolution, much higher in [95, 94]. 
Non-axisymmctric two-dimensional studies, restricted to the equatorial plane of 
the black hole, were discussed in [94], motivated by the non-stationary patterns 
found in Newtonian simulations (see, e.g., [31]). The relativistic computations 
revealed that initially asymptotic uniform flows always accrete on to the hole in 
a stationary way which closely resembles the previous axisymmctric patterns. 

In [219] Papadopoulos and Font presented a procedure which simplifies the 
numerical integration of the general relativistic hydrodynamic equations near 
black holes. This procedure is based on identifying classes of coordinates in 
which the black hole metric is free of coordinate singularities at the horizon (un- 
like the commonly adopted Boyer-Lindquist coordinates), independent of time, 
and admits a spacelike decomposition. With those coordinates the innermost 
radial boundary can be placed inside the horizon, allowing for an unambiguous 
treatment of the entire (exterior) physical domain. In [98, 99] this approach was 
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Figure 9: Relativistic wind accretion onto a rapidly rotating Kerr black hole 
(a = 0.999Ajf, the hole spin is counter-clock wise) in Kerr-Schild coordinates 
(left panel). Isocontours of the logarithm of the density are plotted at the 
final stationary time t — 500Af. Brighter colors (yellow- white) indicate high 
density regions while darker colors (blue) correspond to low density zones. The 
right panel shows how the flow solution looks like when transformed to Boyer- 
Lindquist coordinates. The shock appears here totally wrapped around the 
horizon of the black hole. The box is 12M units long. The simulation employed 
a (r, (/))-grid of 200 x 160 zones. Further details are given in [98]. 

applied to the study of relativistic wind accretion onto rapidly rotating black 
holes. The effects of the black hole spin on the flow morphology were found to 
be confined to the inner regions of the black hole potential well. Within this 
region, the black hole angular momentum drags the flow, wrapping the shock 
structure around. An illustrative example is depicted in Fig. 9. The left panel 
of this figure corresponds to a simulation employing the Kerr-Schild form of 
the Kerr metric, regular at the horizon. The right panel shows how the accre- 
tion pattern would look like were the computation performed using the more 
common Boyer-Lindquist coordinates. The transformation induces a noticeable 
wrapping of the shock around the central hole. The shock would wrap infinitely 
many times before reaching the horizon. As a result, the computation in these 
coordinates would be much more challenging than in Kerr-Schild coordinates. 
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Figure 10: Stills from a movie showing the time evolution of the accre- 
tion/collapse of a quadrupolar shell onto a Schwarzshild black hole. The left 
panel shows isodensity contours and the right panel the associated gravitational 
waveform. The shell, initially centered at r* = 35M, is gradually accreted by 
the black hole, a process which perturbs the black hole and triggers the emission 
of gravitational radiation. After the burst, the remaining evolution shows the 
decay of the black hole quasi- normal mode ringing. By the end of the simulation 
a spherical accretion solution is reached. Further details are given in Rcf. [221]. 
(To see the movie, please go to the electronic version of this review article at 
http : //www . livingreviews . org/Articles/Volume3/2000-2f ont . ) 



4.2.4 Gravitational radiation 

Semi-analytical studies of finite-sized collections of dust, shaped in the form of 
stars or shells, falling isotropically onto a black hole are available in the liter- 
ature [199, 122, 256, 214, 228]. These investigations approximate gravitational 
collapse by a dust shell of mass m falling into a Schwarzschild black hole of mass 
M ^ m. These studies have shown that for a fixed amount of infalling mass 
the gravitational radiation efficiency is reduced compared to the point particle 
limit {AE ^ 0.0104m^/M), never exceeding that of a particle with the same 
mass, the reason being cancellations of the emission from distinct parts of the 
extended object. 

In [221] such conclusions were corroborated with numerical simulations of 
the gravitational radiation emitted during the accretion process of an extended 
object onto a black hole. The first order deviations from the exact black hole 
geometry were approximated using curvature perturbations induced by matter 
sources whose non-linear evolution was integrated using a (non-linear) hydrody- 
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namics code (adopting the conservative formulation of Section 2.1.3 and HRSC 
schemes). All possible types of curvature perturbations are captured in the 
real and imaginary parts of the Weyl tensor scalar (sec, e.g., [59]). In the 
framework of the Ncwman-Penrose formalism the equations for the perturbed 
Weyl tensor components decouple, and when written in the frequency domain 
even separate [285]. Papadopoulos and Font [221] used the limiting case for 
Schwarzschild black holes, i.e., the inhomogeneous Bardeen-Press equation [27]. 
The simulations showed the gradual excitation of the black hole quasi-normal 
mode frequency by sufficiently compact shells. 

An example of these simulations appears in the MPEG movie of Fig. 10. 
This movie shows the time evolution of the shell density (left panel) and the 
associated gravitational waveform during a complete accretion/collapse event. 
The (quadrupolar) shell is parametrized according to p 

with k = 2, po = 10~^ and ro = 35M. Additionally r* denotes a logarithmic 

radial (Schwarzschild) coordinate. The animation shows the gradual collapse 
of the shell on to the black hole. This process triggers the emission of grav- 
itational radiation. In the movie it can be clearly seen how the burst of the 
emission coincides with the most dynamical accretion phase, when the shell 
crosses the peak of the potential and is subsequently captured by the hole. The 
gravitational wave signal coincides with the quasinormal ringing frequency of 
the Schwarzschild black hole, 17M. The existence of an initial burst, separated 
in time from the physical burst, is also noticeable in the movie. It just reflects 
the gravitational radiation content of the initial data (see [221] for a detailed 
explanation) . 

One-dimensional numerical simulations of a self-gravitating perfect fluid ac- 
creting onto a black hole were presented in [223], where the effects of mass 
accretion during the gravitational wave emission from a black hole of growing 
mass were explored. Using the conservative formulation outlined in Section 
2.2.2 and HRSC schemes, Papadopoulos and Font [223] performed the simula- 
tions adopting an ingoing null foliation of a spherically symmetric black hole 
spacetime [222]. Such a foliation penetrates the black hole horizon, allowing 
for an unambiguous numerical treatment of the inner boundary. The essence of 
non-spherical gravitational perturbations was captured by adding to the (char- 
acteristic) Einstein-perfect fluid system, the evolution equation for a minimally 
coupled massless scalar field. The simulations showed the familiar damped- 
oscillatory radiative decay, with both decay rate and frequencies being modu- 
lated by the mass accretion rate. Any appreciable increase in the horizon mass 
during the emission reflects on the instantaneous signal frequency, /, which 
shows a prominent negative branch in the /(/) evolution diagram. The features 
of the frequency evolution pattern reveal key properties of the accretion event, 
such as the total accreted mass and the accretion rate. 

Recently, Zanotti et al. [317] have performed hydrodynamical simulations 
of constant angular momentum thick disks (of typical neutron star densities) 
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orbiting a Schwarzschild black hole. Upon the introduction of perturbations, 
these systems either become unstable to the runaway instability [91] or exhibit 
a regular oscillatory behaviour resulting in a quasi-periodic variation of the 
accretion rate as well as of the mass quadrupole (animations can be visualized 
at the web site listed in Ref. [287]). Zanotti et al. [317] have found that the 
latter is responsible for the emission of intense gravitational radiation whose 
amplitude is comparable or larger than the one expected in stellar core collapse. 
The strength of the gravitational waves emitted and their periodicity are such 
that signal-to-noise ratios ^ 0{1) — 0(10) can be reached for sources at 20 
Mpc or 10 Kpc, respectively, making these new sources of gravitational waves 
potentially detectable. 

4.3 Hydro dynamical evolution of neutron stars 

The numerical investigation of many interesting astrophysical processes involv- 
ing neutron stars, such as the rotational evolution of proto-neutron stars, which 
can be affected by a dynamical bar-mode instability and by the Chandrasekhar- 
Friedman-Schutz instability, or the gravitational radiation from unstable pulsa- 
tion modes or, more importantly, from the catastrophic coalescence and merger 
of neutron star compact binaries, requires the ability of accurate, long-term 
hydrodynamical evolutions employing relativistic gravity. These scenarios are 
receiving increasing attention in recent years [266, 171, 100, 184, 257, 261, 264, 
101, 93, 259]. 

4.3.1 Pulsations of relativistic stars 

The use of relativistic hydrodynamical codes to study the stability properties 
of neutron stars and to compute mode frequencies of oscillations of such ob- 
jects is increasing in recent years (sec, e.g. the Living Reviews article by Ster- 
gioulas [278] and references therein). An obvious advantage of the hydrody- 
namical approach is that it includes, by construction, nonlinear effects, which 
are important in situations where the linearized equations, commonly employed 
in calculations of mode-frequencies of pulsating stars, break down. In addition, 
due to the current status of hydrodynamics codes, the computation of mode fre- 
quencies in rapidly rotating relativistic stars might be easier to achieve through 
nonlinear numerical evolutions than using perturbative computations (see, e.g. 
the results in [93, 259]). 

Hydrodynamical evolutions of polytropic models of spherical neutron stars 
can be used as test-bed computations for multidimensional codes. Representa- 
tive examples are the simulations by [115], with pseudo-spectral methods, and 
by [244] with HRSC schemes. These investigations adopted radial gauge polar 
slicing coordinates in which the general relativistic equations are expressed in a 
simple way which resembles Newtonian hydrodynamics. Gourgoulhon [115] used 
a numerical code to detect, dynamically, the zero value of the fundamental mode 
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of a neutron star against radial oscillations. Romero et al. [244] highlighted the 
accuracy of HRSC schemes by finding, numerically, a change in the stability 
behavior of two slightly different initial neutron star models: for a given EOS, 
a model with mass 1.94532M0 is stable and a model of 1.94518M0 is unstable. 
More recently, in [274] a method based on the non-linear evolution of deviations 
from a background stationary equilibrium star was applied to study nonlinear 
radial oscillations of a neutron star. The accuracy of the approach permitted a 
detailed investigation of nonlinear features such as quadratic and higher order 
mode coupling and nonlinear transfer of energy. 

Axisymmetric pulsations of rotating neutron stars can be excited in several 
scenarios, such as core-collapse, crust and core-quakes and binary mergers, and 
could become detectable either in gravitational waves or high-energy radiation. 
An observational detection of such pulsations would yield valuable information 
about the EOS of relativistic stars. As a first step towards the study of pulsa- 
tions of rapidly-rotating relativistic stars, Font, Stergioulas and Kokkotas [101] 
developed an axisymmetric numerical code which integrates the equations of 
general relativistic hydrodynamics in a fixed background spacetime. The finite- 
difference code is based on a state-of-the-art approximate Riemann solver [74] 
and incorporates different second- and third-order TVD and ENO numerical 
schemes. This code is capable of accurately evolving rapidly rotating stars for 
many rotational periods, even for stars at the mass-shedding limit. The test 
simulations reported in [101] showed that, for non-rotating stars, small ampli- 
tude oscillations have frequencies that agree below the 1% level with linear, 
radial and non-radial, normal mode frequencies in the so-called Cowling ap- 
proximation (i.e. when the evolution of the spacetime variables is neglected). 
Axisymmetric modes of pulsating non rotating stars are computed in [269] both 
in Cowling and fully coupled evolutions. Contrary to the 2+1 approach fol- 
lowed by [101], the code used in Ref. [269] evolves the relativistic stars on null 
spacetime foliations (see Section 2.2.2). 

Until very recently (sec below), the quasi-radial modes of rotating relativis- 
tic stars had been studied only under simplifying assumptions, such as in the 
slow-rotation approximation or in the relativistic Cowling approximation. An 
example of the latter is presented in [92], where using the code developed by 
[101], a comprehensive study of all low-order axisymmetric modes of uniformly 
and rapidly rotating relativistic stars was presented. The frequencies of quasi- 
radial and non-radial modes with spherical harmonic indices £ — 0,1,2 and 
3 were computed through Fourier transforms of the time evolution of selected 
fluid variables. This was done for a sequence of appropriately perturbed sta- 
tionary rotating stars, from the non-rotating limit to the mass-shedding limit. 
The frequencies of the axisymmetric modes are affected significantly by rota- 
tion only when the rotation rate exceeds about 50% of the maximum allowed. 
As expected, at large rotation rates, apparent mode crossings between different 
modes appear. 

In [93], the first mode frequencies of uniformly rotating stars in full gen- 
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eral relativity and rapid rotation were obtained, using the three-dimensional 
code GR_ASTRO described in Section 3.3.2. Such frequencies were computed 
both in fixed spacetimc evolutions (Cowling approximation) and in coupled hy- 
drodynamical and spacctime evolutions. The simulations used a sequence of 
(perturbed) N = 1, K = 100 {G = c = Mq = 1) polytropes of central den- 
sity pc = 1.28 X 10"^, in which the rotation rate varies from zero to 97% 
of the maximum allowed rotational frequency, ^Ix = 0.5363 x 10'' s~^. The 
Cowling runs allowed a comparison with earlier results reported by [92] , obtain- 
ing agreement at the 0.5% level. The fundamental mode-frequencies and their 
first overtones obtained in fully coupled evolutions show a dependence on the 
increased rotation which is similar to the one observed for the corresponding 
frequencies in the Cowling approximation [92]. 

Relativistic hydrodynamical simulations of nonlinear r-modes are presented 
in [279] (see also [158] for Newtonian simulations). The gravitational radiation 
reaction driven instability of the r-modes might have important astrophysi- 
cal implications, provided that the instability does not saturate at low ampli- 
tudes by nonlinear effects or by dissipative mechanisms. Using a version of the 
GR_ASTRD code, Stergioulas and Font [279] found evidence that the maximum 
r-mode amplitude in isentropic stars is of order unity. The dissipative mecha- 
nisms proposed by different authors to limit the mode amplitude include shear 
and bulk viscosity, energy loss to a magnetic field driven by differential rotation, 
shock waves, or the slow leak of the r-mode energy into some short wavelength 
oscillation modes (see [20] and references therein) . The latter mechanism would 
dramatically limit the r-mode amplitude to a small value, much smaller than 
those found in the simulations of [279, 158] (see [278] for a complete list of ref- 
erences on the subject). Energy leak of the r-mode into other fluid modes has 
been recently considered by [117] through Newtonian hydrodynamical simula- 
tions, finding a catastrophic decay of the amplitude only once it has grown to 
a value larger than that reported by [20]. 

The bar-mode instability in differentially rotating stars in general relativity 
has been analyzed by [260] by means of 3+1 hydrodynamical simulations. Using 
the code discussed in Section 3.3.1, Shibata et al. [260] found that the critical 
ratio of rotational kinetic energy to gravitational binding energy for compact 
stars with M/R > 0.1 is ^ 0.24 — 0.25, slightly below the Newtonian value 
~ 0.27 for incompressible Maclaurin spheroids. All unstable stars are found to 
form bars on a dynamical timescale. 

4.3.2 Binary neutron star coalescence 

Many of the current efforts in code development in relativistic astrophysics aim 
at the simulation of the coalescence of compact binaries. Neutron star binaries 
are among the most promising sources of gravitational radiation to be detected 
by the various ground-based interferometers worldwide. The computation of 
the gravitational waveform during the most dynamical phase of the coalescence 
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and plunge depends crucially on hydrodynamical, finite-size eflFects. This phase 
begins once the stars, initially in quasi-equilibrium orbits of gradually smaller 
orbital radius, due to the emission of gravitational waves, reach the so-called 
innermost stable circular orbit. About ~ 10* years after formation of the bi- 
nary system the gravitational wave frequency enters the LIGO/VIRGO high 
frequency band. The final plunge of the two objects takes place on a dynamical 
timescale of a few ms. During the last 15 minutes or so until the stars finally 
merge, the frequency is inside the LIGO/VIRGO sensitivity range. About 16000 
cycles of waveform oscillation can be monitored, while the frequency gradually 
shifts from ~ 10 Hz to ~ 1 kHz. A perturbative treatment of the gravitational 
radiation in the quadrupole approximation is valid as long as M/R <C 1 and 
M/r <Si 1 simultaneously. Ad being the total mass of the binary, R the neutron 
star radius and r the separation of the two stars. As the stars approach each 
other and merge, both inequalities are less valid and fully relativistic hydrody- 
namical calculations become necessary. 

The accurate simulation of a binary neutron star coalescence is, however, one 
of the most challenging tasks in numerical relativity. These scenarios involve 
strong gravitational fields, matter motion with (ultra-) relativistic speeds, and 
relativistic shock waves. The numerical difficulties are exacerbated by the intrin- 
sic multidimensional character of the problem and by the inherent complexities 
in Einstein's theory of gravity, such as coordinate degrees of freedom and the 
possible formation of curvature singularities (e.g. collapse of matter configura- 
tions to black holes). It is thus not surprising that most of the (few) available 
simulations have been attempted in the Newtonian (and post-Newtonian) frame- 
work (see [236] for a review). Many of these studies employ Lagrangian particle 
methods such as SPH, and only a few have considered (less viscous) high-order 
finite difference methods such as PPM [246]. 

Concerning relativistic simulations Wilson's formulation of the hydrody- 
namic equations (see Section 2.1.2) was used in Refs. [303, 304, 172]. Such inves- 
tigations assumed a conformally flat 3-metric, which reduces the (hyperbolic) 
gravitational field equations to a coupled set of elliptic (Poisson-like) equations 
for the lapse function, the shift vector, and the conformal factor. These early 
simulations revealed the unexpected appearance of a "binary-induced collapse 
instability" of the neutron stars, prior to the eventual collapse of the final merged 
object. This effect was reduced, but not eliminated fully, in revised simulations 
[172], after Flanagan [89] pointed out an error in the momentum constraint 
equation as implemented by Wilson and coworkers [303, 304]. A summary of 
this controversy can be found in [236]. Subsequent numerical simulations with 
the full set of Einstein equations (see below) did not find this eS'ect. 

Nakamura and co-workers have been pursueing a programme to simulate 
neutron star binary coalescence in general relativity since the late 1980's (see, 
e.g., [197]). The group developed a three-dimensional code which solves the full 
set of Einstein equations and self- gravitating matter fields [215]. The equations 
are finite-differenced in a uniform Cartesian grid using van Leer's scheme [291] 
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Figure 11: Still from a movie showing the animation of a head-on collision sim- 
ulation of two IAMq neutron stars obtained with a relativistic code [100, 184]. 
The movie shows the evolution of the density and internal energy. The for- 
mation of the black hole in prompt timescales is demonstrated by the sud- 
den appearance of the apparent horizon at t = 0.16 ms {t = 63.194 in code 
units), which is indicated by violet dotted circles representing the trapped pho- 
tons. (To see the movie, please go to the electronic version of this review 
article at http : //www . livingreviews . org/Articles/Volume3/2000-2f ont . ) 
See [3] for download options of higher quality versions of this movie. 

with TVD flux limiters. Shock waves are spread out using a tensor artificial 
viscosity algorithm. The hydrodynamic equations follow Wilson's Eulerian for- 
mulation and the ADM formalism is adopted for the Einstein equations. This 
code has been tested by the study of the gravitational collapse of a rotating 
polytrope to a black hole (comparing to the axisymmetric computation of Stark 
and Piran [276]). Further work to achieve long term stability in simulations of 
neutron star binary coalescence is under way [215]. We note that the hydrody- 
namics part of this code is at the basis of Shibata's code (Section 3.3.1) which 
has successfully been applied to simulate the binary coalescence problem (see 
below) . 

The head-on collision of two neutron stars (a limiting case of a coalescence 
event) was considered by Miller et al. [184], who performed time-dependent rela- 
tivistic simulations using the code described in Section 3.3.2. These simulations 
analyzed whether the collapse of the final object occurs in prompt timescales 
(a few milliseconds) or delayed (after neutrino cooling) timescales (a few sec- 
onds). In [254] it was argued that in a head-on collision event, sufficient thermal 
pressure is generated to support the remnant in quasi-static equilibrium against 
(prompt) collapse prior to slow cooling via neutrino emission (delayed collapse). 
In [184], prompt collapse to a black hole was found in the head-on collision 
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of two IAMq neutron stars modeled by a polytropic EOS with F = 2 and 
K = 1.16 X 10^ cm^/gs^. The stars, initially separated by a proper distance 
of d = 44 km, were boosted towards one another at a speed of GM/d (the 
Newtonian infall velocity). The simulation employed a Cartesian grid of 192'^ 
points. The time evolution of this simulation can be followed in the Quicktime 
movie in Fig. 11. This animation simultaneously shows the rest-mass density 
and the internal energy evolution during the on-axis collision. The formation 
of the black hole in prompt timescales is signalled by the sudden appearance 
of the apparent horizon at i = 0.16 ms (t = 63.194 in code imits). The violet 
dotted circles indicate the trapped photons. The animation also shows a mod- 
erately relativistic shock wave (Lorentz factor ^ 1.2) appearing at i ^ 36 (code 
units; yellow- white colors), which eventually is followed by two opposite moving 
shocks (along the infalling z direction) which propagate along the atmosphere 
surrounding the black hole. 

The most advanced simulations of neutron star coalescence in full general 
relativity are those performed by Shibata and Uryu [257, 264, 265]. Their nu- 
merical code, briefly described in Section 3.3.1, allows the long-term simulation 
of the coalescences of both, irrotational and corotational binaries, from the in- 
nermost stable circular orbit up to the formation and ringdown of the final 
collapsed object (either a black hole or a stable neutron star), includes an ap- 
parent horizon finder, and can extract the gravitational waveforms emitted in 
the collisions. Shibata and Uryu have performed simulations for a large sam- 
ple of parameters of the binary system, such as the compactness of the (equal 
mass) neutron stars (0.12 < M/R < 0.16), the adiabatic index of the F-law 
EOS (1.8 < F < 2.5), and the maximum density, rest mass, gravitational mass 
or total angular momentum. The initial data correspond to quasi-equilibrium 
states, either corotational or irrotational. the latter being more realistic from 
considerations of viscous versus gravitational radiation timescales. These initial 
data are obtained by solving the Einstein constraint equations and the equations 
for the gauge variables luidc^r the assumption of a conformally flat 3-metric and 
the existence of a helicoidal Killing vector (see [265] for a detailed explanation). 
The binaries are chosen at the innermost orbits for which the Lagrange points 
appear at the innca- edge of the neutron stars, and the plunge is induced by 
reducing the initial angular momentum by ^ 2 — 3%. 

The comprehensive parameter space survey carried out by [257, 264, 265] 
shows that the final outcome of the merger depends sensitively on the initial 
compactness of the neutron stars before plunge. Hence, depending on the stiff- 
ness of the EOS, controlled through the value of F, if the total rest mass of the 
system is ~ 1.3 — 1.7 times larger than the maximum rest mass of a spherical 
star in isolation, the end product is a black hole. Otherwise, a marginally-stable 
massive neutron star forms. In the latter case the star is supported against 
self-gravity by rapid differential rotation. The star could eventually collapse 
to a black hole once sufficient angular momentum has dissipated via neutrino 
emission or gravitational radiation. The different outcome of the merger is 
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Figure 12: Snapshots of the density contours in the equatorial plane for a binary 
neutron star coalescence which leads to a rotating black hole (see [265] for the 
characteristics of the model). Vectors indicate the local velocity field (11^,1;^). 
Pt=o denotes the orbital period of the initial configuration. The thick solid circle 
in the final panel indicates the apparent horizon. The figure is taken from [265] 
(used with permission). 



reflected in the gravitational waveforms [265]. Therefore, future detection of 
high-frequency gravitational waves could help to constraint the maximum al- 
lowed mass of neutron stars. In addition, for prompt black hole formation, a 
disk orbiting around the black hole forms, with a mass less than 1% the total 
rest mass of the system. Disk formation during binary neutron star coalescence, 
a fundamental issue for cosmological models of short duration GRBs, is en- 
hanced for unequal mass neutron stars, in which the less massive one is tidally 
disrupted during the evolution (Shibata, private communication). 

A representative example of one of the models simulated by Shibata and 
Uryu is shown in Fig. 12. This figure is taken from Ref. [265]. The compactness 
of each star in isolation is M/R = 0.14 and F = 2.25. Additional properties 
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of the initial model can be found in Table 1 of [265]. The figure shows nine 
snapshots of density isocontours and the velocity field in the equatorial plane 
(z = 0) of the computational domain. At the end of the simulation a black 
hole has formed, as indicated by the thick solid circle in the final snapshot, 
representing the apparent horizon. The formation timescale of the black hole is 
larger the smaller the initial compactness of each star. The snapshots depicted in 
Fig. 12 show that once the stars have merged the object starts oscillating quasi- 
radially before the complete collapse takes place, the lapse function approaching 
zero not monotonically [265]. The collapse toward a black hole sets in after the 
angular momentum of the merged object is dissipated through gravitational 
radiation. Animations of various simulations (including this example) can be 
found at Shibata's website [267]. 

To close this section we mention the work of Duez et al. [76] where, through 
analytic modelling of the inspiral of corotational and irrotational relativistic bi- 
nary neutron stars as a sequence of quasi-equilibrium configurations, the gravi- 
tational wave-train from the last phase (a few hundred cycles) of the inspiral is 
computed. These authors further show a practical procedure to construct the 
entire wave-train through coalescence by matching the late inspiral waveform 
to the one obtained by fully relativistic hydrodynamical simulations as those 
discussed in the above paragraphs [264]. Detailed theoretical waveforms of the 
inspiral and plunge as those reported by [76] arc crucial to enhance the chances 
of experimental detection in conjunction with matched-filtering techniques. 

5 Additional information 

This last section of the review contains technical information concerning re- 
cent developments for the implementation of Riemann solver based numerical 
schemes in general relativistic hydrodynamics. 

5.1 Riemann problems in locally Minkowskian coordinates 

In [230] a procedure to integrate the general relativistic hydrodynamic equations 
(as formulated in Section 2.1.3), taking advantage of the multitude of Riemann 
solvers developed in special relativity, was presented. The approach relies on 
a local change of coordinates in terms of which the spacetime metric is locally 
Minkowskian. This procedure allows, for ID problems, the use of the exact 
solution of the special relativistic Riemann problem [168]. 

Such a coordinate transformation to locally Minkowskian coordinates at each 
numerical interface, assumes that the solution of the Riemann problem is the one 
in special relativity and planar symmetry. This last assumption is equivalent 
to the approach followed in classical fluid dynamics, when using the solution 
of Riemann problems in slab symmetry for problems in cylindrical or spherical 
coordinates, which breaks down near the singular points (e.g., the polar axis in 
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cylindrical coordinates). In analogy to classical fluid dynamics, the numerical 
error depends on the magnitude of the Christoffel symbols, which might be large 
whenever huge gradients or large temporal variations of the gravitational field 
are present. Finer grids and improved time advancing methods will be required 
in those circumstances. 

Following [230] wc illustrate the procedure for computing the second flux 
integral in Eq. (52), which we call T. We begin by expressing the integral 
on a basis with eg = and e| forming an orthonormal basis in the plane 
orthogonal to with ej normal to the surface and eg and tangent to that 
surface. The vectors of this basis verify • = 77.^ with 77^^ the Minkowski 
metric (in the following, caret subscripts will refer to vector components in this 
basis). 

Denoting by Xq the coordinates at the center of the interface at time i, we 
introduce the following locally Minkowskian coordinate system a;" = (a;" — 
Xq), where the matrix M" is given by da = M"ea, calculated at Xq. In this 
system of coordinates the equations of general relativistic hydrodynamics trans- 
form into the equations of special relativistic hydrodynamics, in Cartesian co- 
ordinates, but with non-zero sources, and the flux integral reads 



J V^F^dx°dx'^dx^ = J lF^-^F^]^dx°dx^dx^, (73) 



(the caret symbol representing the numerical flux in Eq. (52) is now removed 
to avoid confusion) with y/—g = 1 + 0{x°'), where we have taken into account 

that, in the coordinates x", E3.1 is described by the equation x^ — ^x^ = 
(with /3* = M?/3'), where the metric elements and a are calculated at Xq. 
Therefore, this surface is not at rest but moves with speed (i^ /a. 

At this point all the theoretical work on special relativistic Riemann solvers 
developed in recent years, can be exploited. The quantity in parenthesis in 
Eq. (73) represents the numerical flux across E^,! , which can, now, be calculated 
by solving the special relativistic Riemann problem defined with the values at 
the two sides of E^,! of two independent thermodynamical variables (namely, 
the rest mass density p and the specific internal energy e) and the components 
of the velocity in the orthonormal spatial basis f* (f* = M^v^). 

Once the Riemann problem has been solved, we can take advantage of the 
self-similar character of the solution of the Riemann problem, which makes it 
constant on the surface E^-i simplifying the calculation of the above integral 
enormously: 



F^-i^F''] / yf^gdx^dx^dx\ (74) 




where the superscript (*) stands for the value on E^-i obtained from the solution 
of the Riemann problem. Notice that the numerical fluxes correspond to the 
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vector fields = {J, T-n, T-ej, T-e2, T-Cg} and linearized Riemann solvers 
provide the numerical fluxes as defined in Eq. (73). Thus the additional relation 

T ■ di = (T ■ Gj) has to be used for the momentum equations. The integral 
in the right hand side of Eq. (74) is the area of the surface T,^! and can be 
expressed in terms of the original coordinates as 

/ sf^gdx^dx^dx^ = [ ^ dx"^ dx^ dx"^ , (75) 

which can be evaluated for a given metric. The interested reader is addressed 
to [230] for details on the testing and calibration of this procedure. 

5.2 Chciracteristic fields in the 3+1 conservative Eulerian 
formulation of Section 2.1.3 

This section collects all information concerning the characteristic structure of 
the general relativistic hydrodynamic equations in the Eulerian formulation of 
Section 2.1.3. As explained in Section 3.1.2 this information is necessary in order 
to implement approximate Riomann solvers in HRSC finite difference schemes. 

We only present the characteristic speeds and fields corresponding to the 
a;-direction. Equivalent expressions for the two other directions can be easily 
obtained by symmetry considerations. The characteristic speeds (eigenvalues) 
of the system are given by: 

Ao = au^-/3^ (triple), (76) 

A± = ^— ^ {v%l - cl) ± c«V(l - ^^2)[7^^(1 - v2c2) - v^v^{l - c2)]} 

- r, ' (77) 

where Cg denotes the local sound speed, which can be obtained from /ic^ — 
X + ^'^) with X = |2 and k = We note that the Minkowskian limit of 
these expressions is recovered properly (see [73]) as well as the Newtonian one 
(Ao = w^, A± = ± Cs). 

A complete set of right-eigenvectors is given by (superscript T denotes trans- 
pose): 



To 



r0,2 = {WVy,h{j^y + 2W'^V^Vy), hijyy + 2W'^ VyVy) , h{j + 2W'^ V ^Vy) , 

Wvy{2hW - l)f , (79) 



ro,3 = {Wv„ h{^^, + 2W\^v,), h{^y, + 2W^VyV,), /i(7^, + 2W\,v,) 



Wv^{2hW-l)f, (80) 
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r± = (1, WC±, hWvy, hWv^, hWA± - if, 
where the following auxiliary quantities are used: 



K — C 



,2 ' 



R = k/p, Cj^ = Va:- V: 



±1 



V3 



Al 



A!t = X±+p\ A = A/a, = f3'/a. 

Finally, a complete set of left-eigenvectors is given by: 



lo,: 



W 



-{h-w, wv"", wyy, Wv\ -wy , 



lo,2 = 



1 



Izzil - V^V") + ^xzVzV'^ 
-lvz{l - VxV^) - IxzVyV'' 

^zz'^y ~t~ ^yz'^z 



lo,3 = 



1 



lyyVz + IzyVy 

V'^i-fyyVz-JzyVy) 

-7zj/(l - VxV'') - Ixy'^zV'- 

-1yy{\ - VxV'') ^ -ixyVyV"" 

—'^yyVz + IzyVy 
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hwvu + 

r^^(l - KAl) + (2/C - l)Vl{W^v^^ - T^^v^) 
T^y{l - 1CA%) + (2/C - l)Vl{W^vy^ - T^yv^) 
r,,(l - ICAl) + (2/C - l)Vl{W^v^S, - T^,v-) 

L (1 - /c) [-71;^ + VI {w^^ - r^^)] - KWvu 

where the following relations and aiixiliary quantities have been used: 



4 



(88) 



l-Al = v^V: 



(89) 



{CI - CI) + {Aiyi - A%VV) = 0, 



(90) 



A=h^W{)C-l){C'l-Cl)^ ^ = r,,-'yv^v'' 



(91) 



(92) 



^xx lyyizz 'lyz' ^xy 'Jyzlxz Txj/Tzz) ^xz Ixylyz ^xzlyy: (^3) 



^xx'^x 
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(94) 



These two sets of eigenfields reduce to the corresponding ones in the Minkowskian 
limit [73]. 
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